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ABSTRACT 


The  Third  Quarterly  Technical  Report  for  the  LASA  Experimental  Signal 
Processing  System  identifies  the  effort  expended  to  provide  the  hardware  and 
software  in  support  of  research  and  development  directed  toward  the  study  of 
seismic  signal  processing.  It  also  delineates  work  tasks  planned  for  the  next 
quarter.  This  document  presents  detailed  information  related  to  machine  con¬ 
figurations,  time  delay  correlation,  event  location  accuracy,  optimum  process¬ 
ing,  fast  Fourier  transform,  array  design,  steering  delay  library,  magnitude 
estimation,  and  travel  time  characterization. 
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Section  1 


INTRODUCTION 


The  work  reported  herein  was  performed  under  the  "LASA  Experimental 
Signal  Processing  Systemn  Contract  Number  F19628-67-C-0198,  and  is  a  con¬ 
tinuation  find  extension  of  the  "Large  Aperture  Seismic  Array  Signal  Processing 
Study, "1  Contract  Number  SD~296,and  the  nLarge  Aperture  Seismic  Array 
Signal  Processing  Communications  and  Simulation  Study, "2,3, 4  Contract  Number 
A  FI  9(62  8) -594 8. 

The  purpose  of  this  contract  is  to  design,  develop,  and  implement  the  LASA 
Experimental  Signal  Processing  System  (ESPS),  including  the  hardware  and  soft¬ 
ware,  to  provide  an  experimental  capability  to: 

a.  Evaluate  performance  of  the  system  in  accordance  with  system 
requirements 

b.  Demonstrate  the  capability  to  meet  basic  LASA  signal  processing 
objectives 

c.  Conduct  research  to  develop  means  of  improving  and  extending 
the  capability  of  the  system 

d.  Perform  seismic  and  signal  processing  experiments  of  interest. 

5  0 

During  the  first'  and  second  quarters,  effort  was  primarily  directed 
toward  establishing  the  Seismic  Array  Analysis  Center  (SAAC),  Washington,  D.C., 
installation  and  operation  of  the  Detection  and  Event  Processors  as  System/360 
Model  G  and  II  general  purpose  digital  computers,  respectively,  and  the  con¬ 
tinuation  of  system  studies  in  seismic  signal  processing. 

Work  in  the  present  quarter  has  focused  on  such  signal  processing  studies 
as  time  delay  correlation,  event  location  accuracy,  optimum  processing,  fixed- 
point  Fourier  transform,  and  array  design  programming  for  system  development. 
In  addition,  activity  in  event  magnitude  estimating,  time  delay  library  formula¬ 
tion,  and  travel  time  characterization  supporting  event  processing  has  been 
initiated. 
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Seetion  2 


RESULTS  AND  SUMMARY  OF  WORK  ACCOMPLISHMENTS 


The  following  subsections  correspond  to  the  functional  disciplines  of  opera¬ 
tions,  systems,  and  programming. 


2.1  OPERATIONS 

The  principal  effort  during  this  quarter  was  expended  on  the  maehine  sys¬ 
tem,  Experimental  Display  and  computer  operations. 

Installation  and  checkout  of  the  SAAC  hardware  configuration  was  completed. 
The  current  and  planned  configurations  are  described  in  Appendix  I.  The  IBM 
2250-1  CRT  Display  was  installed,  and  debug  of  the  software  support  packages 
was  initiated. 

Modification  and  integration  of  the  Experimental  Display  was  completed  for 
both  the  off-line  and  on-line  modes.  The  strip  chart  recorder  interface, 
reported  previously  has  been  deferred  due  to  priority  work  in  other  areas. 


2.2  SYSTEMS 

During  this  quarter,  effort  pertinent  to  overall  process  improvement  and 
event  processing  structure  has  been  undertaken. 

As  reported  in  Appendix  II- 1,  a  method  of  automatically  determining  time 
delays  from  the  event  data  has  been  formulated.  Preliminary  testing  indicates 
process  stability  and  that  the  resulting  set  of  time  delays  is  consistent  even  when 
the  signal-to-noise  ratio  at  the  seismometer  level  is  near  unity. 

The  effects  of  time  delay  and  system  calibration  errors  on  large  array 
event  location  aecuraey  are  presented  in  Appendix  II-2.  It  is  shown  that  although 
the  processed  output  record  signal-to-noise  ratio  and  system  calibration  errors 
contribute  to  location  error,  the  significant  location  errors  are  due  to  array 
aperture  and  incorrect  time  delays. 
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A  study  to  determine  the  effectiveness  of  fidelity  optimum  processing  to 
cJiscriminate  against  one  of  two  simultaneously  occurring  events  has  been  under¬ 
taken.  The  results  presented  in  Appendix  II- 3  relate  the  effect  of  signal-to-noise 
ratio  and  array  geometry  on  processing  performance. 

Current  interest  in  the  Fast  Fourier  Transform  algorithm  has  prompted 
investigation  of  its  fixed  arithmetic  execution  precision  to  determine  the  appli¬ 
cability  of  microprogramming.  It  is  shown  in  Appendix  II-4  that  the  errors  in 
fixed  point  execution  can  be  reasonably  bounded  by  scaling  the  process  when  an 
overflow  occurs  and  that  microcoded  execution  speed  factors  are  about  ten  over 
standard  assembly  language  programming. 

A  method  of  array  geometry  optimization  has  been  established  and  tested 
and  is  described  in  Appendix  II-5.  Since  the  technique  employs  interelement 
distance  dependent  correlation  coefficient  data  as  input,  its  principal  advan¬ 
tage  appears  to  be  comparative  in  nature. 

In  the  phase  delay  activity,  effort  was  concentrated  on  defining  process 
structure,  automatic  acquisition,  anomaly  characterization,  and  subarray  beam 
assignment. 

With  respect  to  process  structure  and  automatic  acquisition,  the  chief 
accomplishment  has  been  the  development  of  a  library  organization  for  handling 
the  phase  delays  for  both  detection  and  event  processing.  Progress  on  this 
effort  is  described  in  Appendix  III- 1. 

A  procedure  based  on  the  use  of  sets  of  correction  factors  where  each  phase 
delay  set  is  valid  over  a  specified  geographic  area  called  a  region  has  been 
developed.  The  procedure  also  involves  the  use  of  interpolation  techniques  to 
calculate  corrections  relative  to  areas  between  regions.  Appendix  III— 2  describes 
the  regions  and  their  associated  correction  factors  as  well  as  a  test  which  was 
performed  to  evaluate  the  effectiveness  of  the  region  groupings.  The  results  of 
the  test  indicate  that  regional  corrections  improve  system  performance. 

The  subarray  beam  assignment  investigation  addresses  the  method  of  cover¬ 
ing  an  arbitrary  location  in  inverse  velocity  space  by  a  judicious  choice  of  sub¬ 
array  beams.  During  this  quarter,  formulas  were  developed  for  calculating  sub¬ 
array  steering  losses  relative  to  a  specified  configuration  of  preformed  beams 
and  the  associated  LASA  steering  loss  resulting  from  the  use  of  such  preformed 
beams.  In  addition,  a  rule  was  specified  for  computing  the  beam  from  each  sub¬ 
array  which  is  to  be  used  in  forming  the  LASA  beam,  when  given  a  desired  LASA 
beam  location.  Details  regarding  the  subarray  assignment  investigation  may  be 
found  in  Appendix  III- 3. 

In  Appendix  IV- 1,  the  empirical  relationships  between  the  classical  event 
magnitude  measurement  and  the  kinetic  energy  "density1'  method  of  magnitude 
estimation  are  given. 
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A  significant  portion  of  event  processing  assumes  the  use  of  travel  time 
data  to  sort  wave  arrivals  into  consistent  arrival  families  and  to  use  this  data 
in  estimating  event  location.  The  travel  time  information  in  its  original  tabular 
form  is  bulky  and  would  require  interrelation  with  epicentral  range  and  haz¬ 
ardous  differentiation  of  discrete  data.  In  Appendix  IV-2,  orthogonal  polynomial 
fitting  in  the  least  squares  error  sense  was  applied  to  a  portion  of  the  Jeffreys- 
Bullen  Seismological  Tables  to  explore  the  existence  of  explicit  functions  relating 
body  wave  travel  time  and  horizontal  inverse  phase  velocity  to  epicentral  event 
range  at  a  known  depth. 


2.3  PROGRAMMING 

Programming  effort  during  this  quarter  was  concentrated  in  Detection  Pro¬ 
cessing,  Supervisor/Monitor,  Plotting  Routines,  Display  Programs,  and  Experi¬ 
mental  Programs.  Accomplishments  in  each  of  these  areas  will  be  presented  in 
the  following  subsections. 


2.3.1  Detection  Processing 

At  the  completion  of  this  quarter's  work,  three  operational  versions  of  the 
Detection  Controller  were  in  existence.  Two  of  these  were  oriented  towards  a 
real-time  operating  environment,  while  the  third  was  tailored  to  support  analysis 
activity. 

The  difference  between  the  two  real-time  versions  lies  in  the  area  of  input/ 
output  control.  One  operates  using  the  standard  IBM  System /.3(30  Disk  Oper¬ 
ating  System  (DOS),  while  the  other  operates  using  a  specially  written  input/ 
output  controller.  Some  experiments  were  made  to  determine  if  there  was  any 
timing  advantage  of  one  method  or  the  other.  The  latter  was  faster  by  about 
0.5  percent  and  appears  to  have  no  significant  effect  on  the  total  number  of  array 
beams  that  can  be  performed  in  a  given  time  interval. 

The  data  analysis  version  of  the  Detection  Controller  uses  standard  DOS 
input/output,  but  it  is  intended  for  operation  under  slightly  different  conditions. 

In  particular,  the  signal  rectification  and  integration  process  has  been  set  at 
5  Hz,  rather  than  1.67  Hz,  and  the  threshold  detection  and  noise  integration 
functions  have  been  removed.  The  capability  to  produce  intermediate  tapes  of 
subarray  beams  (unfiltered  and  filtered),  array  beams,  and  short  time  average 
values  has  been  included  in  this  version.  The  primary  use  of  this  program, 
thus  far,  has  been  to  produce  input  tapes  for  both  beam  pattern  correlation  and 
experimental  display  work. 

All  three  of  these  programs  have  the  same  basic  structure  and  design 
philosophy.  All  accept  the  standard  SAAC  Edit  Format,  as  input,  although  they 
do  not  presently  interpret  the  array  status  data.  A  unique  feature  of  these  pro¬ 
grams  is  a  storage  allocation  routine  which  allows  maximum  use  of  available 
core  storage  for  each  run  of  the  Detection  Controller.  This  has  enabled  us  to 
process  as  many  as  1024  LASA  beams  in  one  run  of  the  program. 
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An  additional  pi'Ogram  capability  includes  the  dispersed  subarray  beam 
capability  (see  Appendix  III- 3).  An  input  card  supplied  with  the  steering 
delays  for  each  array  beam  indicates  to  the  Detection  Controller  which  of  the 
subarray  beams  are  to  be  used  for  array  beamforming.  On  the  basis  of  this 
information,  the  delays  are  properly  indexed  for  use  by  the  beamforming 
microprogram. 


2.3.2  Supervisor/Monitor 

Systems  generation  and  maintenance  activity  was  accented  upon  availability 
of  the  System/360  2040H  early  in  the  quarter.  A  multiprogramming  Disk  Opera¬ 
ting  System  was  provided.  Various  versions  were  generated  as  configuration 
changes  took  place.  At  this  point,  the  system  contained  no  alterations,  but  it  was 
generated  with  a  greatly  expended  supervisor  area  in  anticipation  of  the  future 
core  memory  requirement.  Special  attention  was  given  to  output  spooling  tech¬ 
niques  which  would  allow  reasonably  efficient  program  development  using  a 
multisystem  including  just  one  'line  printer  and  one  card  punch. 


2. 3.2.1  Detection  Subsystem 

A  test  which  measured  the  overhead  incurred  by  using  a  standard  Disk 
Operating  System  supervisor  in  a  detection  processing  function  was  conducted. 
The  results  underscored  the  advantages  in  choosing  this  system  as  a  super¬ 
structure  around  which  the  Detection  Monitor  could  be  built. 

Work  in  this  quarter  also  included  detail  design,  coding,  and  simulated 
testing  of  the  first  Detection  Monitor  package.  This  package  embraced  capa¬ 
bilities  originally  planned  as  two  models.  First,  the  Disk  Operating  System  was 
modified  to  enable  multiprogramming  operation  in  a  machine  without  storage 
protection  so  that  the  data  acquisition  function  could  operate  asynchronously 
with  the  detection  processing  function,  while  having  common  access  to  certain 
portions  of  core  memory.  The  monitor  was  designed  so  that  inputs  could  be 
received  in  real  time  from  any  source,  with  the  detection  processing  programs 
relieved  of  device-dependent  considerations.  The  same  design  also  enabled  the 
setting  of  mode  indicators  from  external  sources  to  allow  dynamic  entry  by  the 
detection  processing  program  into  various  outputting  modes  for  display  purposes. 

The  design  approach  is  notable  in  that  the  Detection  Controller  may  be  inte¬ 
grated  with  the  Detection  Monitor  into  a  working  system  merely  by  assembling 
macros  with  the  program.  The  code  generated  by  the  macros  represents  the 
Detection  Monitor  and  becomes  a  separate  program  only  at  the  time  the  Detec¬ 
tion  Subsystem  is  loaded  and  initialized  for  execution.  This  was  tested  by  using 
a  simulated  Detection  Controller  and  a  program  which  was  executed  in  the  other 
processor  to  play  the  role  of  the  real-time  interface  function. 
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As  the  quarter  closed,  work  began  on  a  third  Detection  Monitor  model. 
DEMON  HI  represents  a  structural  rather  than  a  functional  change  over  pre¬ 
vious  versions.  The  Detection  Subsystem  programs  in  DEMON  III  will  be  entirely 
core-resident  at  all  times.  In  previous  versions,  programs  were  loaded  from 
disk  memory  to  perform  infrequent  functions  such  as  handling  I/O  device  errors 
and  operator  communication.  The  delays  introduced  by  such  an  arrangement 
cannot  be  tolerated  in  a  real-time  system.  Therefore  an  appropriate  subset 
of  the  Disk  Operating  System  transient  functions  will  be  made  a  part  of  the 
permanently  resident  nucleus. 


2. 3. 2. 2  Event  Subsystem 

Significant  progress  was  made  toward  providing  operating  sytems  tailored 
specifically  to  SAAC  needs.  These  systems  take  the  form  of  modifications  to 
the  Disk  Operating  System.  The  manner  in  which  the  modifications  are  provided 
allows  the  generation  of  a  standard  DOS,  an  SAAC  tailored  model,  or  anything 
in  between,  by  specifying  parametric  values  associated  with  one  set  of  macro 
source  statements.  That  is,  only  one  system  source  copy  must  be  maintained. 

The  first  step  in  providing  suitable  systems  was  to  include  the  support  of 
devices  in  the  SAAC  configuration  which  are  not  supported  in  DOS.  This  fulfills 
our  needs  for  program  debugging  at  SAAC  and  serves  as  the  basis  upon  which 
other  monitors  may  be  developed. 

These  changes  to  the  supervisor  were  implemented  and  tested  during  this 
quarter,  and  the  package  was  provided  for  daily  general-purpose  operation  in 
place  of  the  standard  system.  It  includes  input/output  macros  and  error  cheeking 
for  the  following'  devices: 

1G27  Plotter 
Channel-to-Channel 
2250  Model  I  Display 

Experimental  Display/Strip  Chart  Recorder 

Also,  an  effort  was  made  to  provide  a  higher  level  of  2250  Display  support 
by  moving  a  capability  from  the  IBM/3G0  Operating  System  into  the  SAAC  system. 

A  fundamental  requirement  for  the  monitor  needed  to  support  development 
of  event  processing  programs  and  to  support  the  ultimate  automatic  operation  is 
the  capability  to  schedule  and  execute  tasks  automatically,  as  previously  out¬ 
lined.  The  capability  to  dynamically  schedule  tasks  into  the  background  job 
queue  has  been  designed  and  coded  as  an  extensive  set  of  modifications  to  the 
Disk  Operating  System.  This  scheduling  program  must  also  coordinate  the 
spooling  of  outputs  generated  by  two  programs  executing  simultaneously  and 
intended  for  the  same  destination  (e.g. ,  printing).  These  functions  were  under¬ 
going  test  as  the  quarter  ended.  An  extension  of  this  will  allow  a  list  of  jobs 
to  be  created  and  managed  dynamically  for  the  foreground  partition  as  well  as 
the  background  partition.  This  work  is  in  the  design  phase. 
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Other  support  includes  the  specification  for  a  possible  remote  terminal 
operation  which  can  be  built  into  the  monitor  discussed  above. 


2.3.3  Plotting  Routines 

The  SAAC  machine  configuration  contains  an  on-line  IBM  1627  plotter. 

The  step  size  of  this  particular  digital  plotter  is  0.01  inch.  The  pen  can  be 
moved  in  three  directions:  horizontally,  vertically  and  diagonally.  Moves 
which  are  not  in  one  of  these  basic  directions  are  accomplished  by  a  combina¬ 
tion  of  the  basic  movements  in  such  a  way  that  the  approximation  is  always 
within  one-half  the  step  size. 

To  facilitate  the  use  of  the  plotter,  several  subroutines  were  written  for 
use  with  FORTRAN  programs.  Instructions  for  the  plotter  are  written  on  mag¬ 
netic  tape  by  these  routines  and  later  sent  to  the  plotter  by  a  utility  program 
which  may  be  run  simultaneously  with  other  programs. 

Three  basic  routines  are  available  for  line  generation  and  annotation  of 
plots.  Their  FORTRAN  names  are  PLOT,  SYMBOL  ,  and  NUMBER. 

The  PLOT  subroutine  calculates  the  necessary  plotter  instructions  to  move 
from  one  coordinate  to  another  with  the  pen  raised  or  lowered.  The  coordinates 
are  given  in  inches  of  deflection  from  a  present  origin.  The  algorithm  to  gener¬ 
ate  lines  which  are  not  in  one  of  the  basic  directions  is  an  integral  part  of  the 
routine  as  is  the  handling  of  all  output  instructions. 

The  SYMBOL  subroutine  accomplishes  the  translation  of  alphameric  infor¬ 
mation  into  sets  of  pen  movements.  Titles  and  other  annotation,  as  well  as 
symbols  drawn  at  specific  data  points,  can  be  drawn  from  lists  of  letters  on 
special  codes. 

Numeric  data  is  translated  to  its  alphameric  equivalent,  with  the  appropriate 
number  of  decimal  places,  by  the  NUMBER  subroutine  and  then  sent  to  the 
SYMBOL  routine. 


2.3.4  Display  Program 

Design  and  implementation  of  an  initial  program  capability  for  both  the 
Experimental  Display  and  the  IBM  2250  Visual  Display  unit  were  accomplished 
during  this  quarter. 

The  2250  display  program  accepts  as  input  the  seismometer  level  edited 
tapes  and  displays  up  to  eight  channels  on  the  display.  The  channels  are 
selectable  through  the  2250fs  console  typewriter.  Scaling  is  controlled  through 
the  functional  keyboard  of  the  2250.  Appropriate  time  and  channel  identification 
information  is  displayed  in  addition  to  the  seismometer  waveforms. 
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The  Experiment  Display  program  accepts  as  input  the  rectified  and  inte¬ 
grated  beam  tapes  produced  by  the  Detection  Controller  and  displays  the  traces 
as  a  32x32  matrix  on  the  cathode  ray  tube  experimental  display.  The  beams 
can  be  displayed  in  a  forward  or  reverse  direction,  in  real  time  or  at  a  fast  scan¬ 
ning  rate,  in  a  single  step  mode  in  which  the  display  is  updated  whenever  an 
interrupt  button  is  pushed,  or  in  a  stationary  time  state.  The  scaling  and  pedes¬ 
tal  values  which  dynamically  control  the  presentation  are  continuously  adjustable 
from  the  Experimental  Display  Unit. 


2.3.5  Experimental  Programs 

The  implementation  of  two  experimental  programs— Time  Delay  Correlation 
and  Beam  Pattern  Correlation— began  during  this  quarter,  and  an  initial  capability 
was  developed  for  both  areas.  Continued  modification  of  these  programs  is  antici¬ 
pated  as  experimental  results  evolve.  They  are  discussed  briefly  in  the  following 
paragraphs. 

The  Time  Delay  Correlation  program  was  written  to  provide  an  automatic 
means  for  determining  subarray  beam  steering  delays.  Using  an  iterative  pro¬ 
cess,  the  program  performs  a  cross- correlation  between  an  array  beam  and  each 
of  the  subarray  beams  or  sensors  in  turn,  adjusting  each  subarray  beam  delay 
based  upon  the  incidence  of  the  maximum  cross-correlation.  A  more  detailed 
description  will  be  found  in  Appendix  II- 1  of  this  report. 

A  Beam  Pattern  Correlation  program  was  developed  to  scan  the  Detection 
Processor  Short  Time  Average  (STA)  output  for  seismic  events.  The  program 
accepts  a  reference  pattern  as  input  along  with  the  STA  data  and  generates  a 
display  tape  as  output  as  well  as  the  printed  output.  The  correlation  is  performed 
over  each  STA  time  sample  of  1024  Beams  to  generate  a  set  of  correlation  indexes. 
These  indexes  represent  the  correlations  for  each  position  of  the  reference  pat¬ 
tern  in  the  beam  field.  These  values  are  then  adjusted  and  placed  on  the  display 
tape  in  such  a  way  that  the  positions  with  the  largest  correlations  will  be  shown  as 
high  intensity  areas  on  the  display.  Each  input  STA  time  sample  generates  one 
display  time  sample.  This  work  is  currently  in  progress. 


8 


Section  3 


PLANS 


The  following  activities  will  be  initiated  during  the  next  quarter.  Discussions 
correspond  to  the  functional  disciplines.  Each  task  is  presented  in  sufficient 
perspective  to  identify  its  interrelation  and  contributions  to  the  project. 


3.1  OPERATIONS 

Activity  for  the  next  reporting  period  centers  upon  completing  plans  for  SAAC 
machine  expansion  to  three  computers,  as  indicated  in  Appendix  1-1. 

Plans  call  for  machine  utilization  to  increase  with  a  partial  second  shift  opera¬ 
tion  scheduled  throughout  the  period.  In  addition  to  program  verification  and  debug 
efforts  continuing  at  an  accelerated  level,  two  other  activities  will  utilize  major 
blocks  of  computer  time.  These  activities  include  the  generation  of  a  modest  library 
of  edited  LASA  event  tapes  and  an  expanded  effort  in  data  analysis  of  these  events. 

The  availability  of  applications  programs  will  allow  investigations  of  the  means 
of  utilization,  limitations,  and  capabilities  of  both  the  IBM  2250  CRT  Display  and 
the  Experimental  Beam  Display.  The  evaluation  of  these  displays  is  expected  to 
have  a  major  effect  on  the  definition  of  requirements  for  system  maintenance  and 
operations  consoles. 


3.2  SYSTEMS 

The  time  delay  correlation  program  reported  in  Appendix  II- 1  appears  to  war¬ 
rant  extension.  Emphasis  will  be  placed  on  experimentation  for  its  intended  pur¬ 
pose  and  upon  applying  the  technique  to  event  location.  In  the  area  of  computa¬ 
tional  signal  processing  techniques,  null  steering  development  and  spatial  beam 
detection  will  be  investigated  as  potential  system  functional  components. 

With  respect  to  the  time  delay  library  process  structure,  the  chief  tasks 
ahead  are  the  following: 

a.  Specify,  in  detail,  a  maintenance  procedure  and  schedule  for  checking 
the  library  of  phase  delays 
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b.  Specify  the  conditions  under  which  the  correlation  program  will  be 
used  for  updating,  correcting,  and  enlarging  the  library 

c.  Establish  criteria  for  identifying  events  to  be  used  in  library 
updating. 

Programming  of  the  sector  interpolation  formulas  will  be  undertaken.  It  is 
hoped  that  the  effectiveness  of  these  formulas  in  obtaining  accurate  steering 
delays  will  be  tested  by  applying  them  to  several  events  located  outside  the  35 
regions. 

The  magnitude-kinetic  energy  "density11  relationship  will  be  used  to  find 
event  magnitudes  for  a  number  of  events  using  unfiltered  and  filtered  LASA 
beams  and  various  window  sizes.  Comparisons  will  be  made  with  event  magni¬ 
tudes  calculated  by  the  classical  method.  The  differences  in  magnitude  esti¬ 
mates  of  each  LASA  subarray  to  that  obtained  on  an  array  basis  should  be 
explored  for  different  values  of  magnitude. 

An  attempt  will  be  made  to  find  polynomials  which  fit  body  wave  travel  time 
and  horizontal  inverse  phase  velocity  as  functions  of  both  epicentral  event  range 
and  event  depth. 

In  addition,  it  appears  desirable  to  investigate  certain  system  instrumentation 
configurations  to  establish  possible  interim  on-line  operating  objectives  and 
requirements. 


3.3  PROGRAMMING 

The  planned  effort  includes  activity  in  the  Detection  and  Event  Processor 
areas,  as  well  as  in  the  continued  development  of  experimental  and  data  analysis 
programs. 

Primary  emphasis  in  the  Detection  Processor  area  will  be  on  extending  the 
existing  capability  of  the  real-time  version  of  the  Detection  Controller  to  operate 
with  varying  array  configurations.  This  extension  will  provide  the  mechanism 
to  evaluate  storage  and  timing  considerations  as  well  as  detection  capability  for 
different  magnitude  events.  In  the  Detection  Monitor  area,  design  and  imple¬ 
mentation  will  proceed  on  DEMON  III,  the  core  resident  version.  A  version  of 
the  Detection  Controller  will  be  merged  with  the  first  Detection  Monitor  Model 
and  tested  by  simulation. 

In  the  Event  Processor  area,  the  present  supervisor  will  be  extended  to 
allow  the  experimental  display  to  function  more  efficiently  with  the  processor. 
Currently,  when  the  experimental  display  program  is  executing,  the  Event 
Processor  must  be  dedicated  to  this  function,  thereby  preventing  the  normal 
multiprogramming  mode  of  operation. 
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In  the  experimental  and  data  analysis  programming  areas,  efforts  will  be 
continued  in  developing  and  extending  the  experimental  and  IBM  2250  display 
capability,  time  delay  correlation,  noise  correlation,  and  beam  pattern  correla¬ 
tion  areas. 
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Appendix  I 


SAAC  COMPUTER  EQUIPMENT 


Figure  1  depicts  the  configuration  of  SAAC  machines  at  the  end  ot  this 
reporting  period.  Changes  made  to  the  configuration  sinee  the  last  report  fall 
into  two  categories: 

a.  Additions  made  to  implement  the  planned  growth  of  the  laboratory. 
Included  in  this  category  are  the  addition  of  the  IBM  2250-1  Display 
Unit,  the  IBM  2740  Remote  Printer  Keyboard,  the  IBM  029  Inter¬ 
preting  Keypunch,  and  the  expanded  capability  Experimental 
Display. 

b.  Additions  and/or  reconfigurations  undertaken  to  better  adapt  the 
system  to  its  evolving  utilization.  Because  of  the  heavy  emphasis 
on  program  evaluation  and  debug,  independent  operation  of  each 
System/360  Model  40  became  increasingly  desirable.  As  a  result 
a  2841-1  Disk  Control  Unit  and  2311-1  Disk  Unit  were  added  to 

the  Detection  Processor  so  that  each  processor  could  operate  under 
an  independent  Disk  Operating  System. 

Integration  of  the  Experimental  Display  with  the  Event  Processor  has 
released  a  2403-2  Tape  Drive  and  Control  Unit.  These  were  connected  to  the 
Detection  Processor  Channel  No.  2  to  provide  a  tape  dump  for  print-out.  The 
resulting  tape  is  printed  by  the  Event  Processor  on  the  1403-2  line  printer 
operating  in  a  foreground  mode  and  overlapped  with  normal  Event  Processor 
operation. 

Figure  2  shows  the  machine  system  developments  planned  during  the  next 
reporting  period.  The  addition  of  the  1403-2  Line  Printer  to  the  40G  will 
eliminate  the  need  for  the  2403  Tape  Control  and  Tape  Unit  and  allow  it  to  be 
deleted  from  the  system.  The  replacement  of  the  2520  by  the  2540  Reader- Punch 
will  improve  the  efficiency  of  card  operations. 

Figure  3  illustrates  the  planned  Auxiliary  System  configuration. 
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Figu  e  1  Co  rent  Eqoipmen*  Configuration 
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Figu  e  2  Planned  Equipment  Configuration 
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Figure  3.  Auxiliary  Equipment  Configuration 


Appendix  II 


SIGNAL  PROCESSING 


In  this  appendix,  the  technique  of  automatic  time  delay  generation  and  the  esti¬ 
mated  event  location  accuracy  of  large  arrays  are  presented.  The  results  of  a 
study  on  the  use  of  fidelity  optimum  processing  against  one  of  two  simultaneously 
occurring  events  appear  promising  and  warrant  both  study  continuation  and  initial 
experimentation.  The  accuracy  of  fixed  point  arithmetic  execution  and  the  esti¬ 
mated  execution  speed  improvement  obtained  by  microprogramming  the  Fast 
Fourier  Transform  is  discussed.  The  results  of  studies  in  array  optimization  are 
also  included. 


II.  1  TIME  DELAY  CORRELATION 


II.  1.1  Introduction 

The  time  delay  correlation  procedure  described  herein  was  suggested  as  a 
topic  of  investigation  because  it  was  felt  that  the  large  volume  of  time  delay  data 
anticipated  during  system  operations  could  not  be  efficiently  and  precisely  handled 
manually.  The  results  of  this  study  and  especially  its  implications  on  event  pro¬ 
cessing  warrant  further  priority  considerations.  In  the  present  section,  the  method 
and  some  experimental  results  are  given  and  the  effects  of  incorporating  wavefront 
analysis  are  also  discussed.  The  potential  effects  of  the  findings  on  system  per¬ 
formance  and  recommendations  are  given  in  the  conclusions. 


II. 1.2  Procedure 

Depicted  in  Figure  4  is  the  formulation  of  the  time  delay  correlation  procedure. 
It  employs  plane  wave  delays  with  regional  corrections  (see  Appendix  III- 2)  as 

initial  values  for  the  time  delays  (Tj_,  t 2, - ,  rn).  The  Laplace  notation  is  used  to 

signify  the  delay  process  from  which  the  beam  [B(t) ]  is  formed.  Each  input  pc  1  (t) , 

X2(t), - ,X  (t)]  is  correlated  with  the  beam.  The  time  (fn)  at  which  the  peak  value 

of  the  correlation  function  [Rn(r)]  occurs  is  selected  and  used  to  estimate  the  cor¬ 
rection  to  the  initial  time  delay  values.  The  quantity  a  has  been  introduced  to 
provide  a  process  stability  control  parameter  (i.e.  when  a  1,  full  value  correc¬ 
tions  are  applied).  In  the  process  programming,  operational  coding  notation  and 
floating  point  arithmetic  have  been  employed. 
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Figure  4.  Time  Delay  Correlation 
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Once  initialized,  the  process  is  allowed  to  iterate  K  times.  As  soon  as 
fn(K^l)  is  not  significantly  different  from  rn(K),  it  is  conjectured  that  the  beam 
power  is  near  maximum  and  the  time  delays  are  optimal. 

For  the  limited  testing  performed  to  date,  the  procedure  has  been  found  to  be 
susceptible  to  false  nulls.  As  expected,  it  has  been  observed  that  the  correlation 
function  is  sinusoidal  and  that  suboptimal  convergence  occasionally  occurs  at  an 
integral  number  of  periods  away  from  maximum.  The  system’s  narrow  band 
response  can  be  employed  to  remove  this  ambiguity  by  using  wave  front  analysis 
techniques  to  test  the  observed  nulls  for  credibility.  Since  the  false  nulls  are  large 
with  respect  to  wave  front  corrugations,  the  procedure  is  expected  to  yield  satis¬ 
factory  results. 


II. 1.3  Experimental  Results 

To  demonstrate  technique,  a  limited  sample  of  certain  experimental  results 
is  given.  Although  use  of  subarray  beams  would  no  doubt  yield  better  results, 
single  seismometers  were  employed  for  computational  convenience.  Depicted  in 
the  upper  portion  of  Figure  5  is  a  well-behaved  element  trace  (a)  and  its  beam 
correlation  function  (b) .  A  poor  element  signal  (c)  and  its  beam  correlation  func¬ 
tion  (d)  are  shown  by  the  lower  pair.  Note  that  the  sensitivity  of  (b)  is  much 
greater  than  (d)  and,  hence,  for  this  data,  the  confidence  of  estimating  the  time 
delay  from  (b)  is  significantly  greater  than  that  for  (d). 

The  effect  of  signal- to- noise  ratio  was  also  investigated  by  applying  the  cor¬ 
relation  procedure  to  scaled  Longshot  data  for  which  the  average  input  signal-to- 
noise  ratio  was  significantly  reduced.  The  average  of  the  magnitude  of  the  differences 
between  the  correlation  time  delays  and  the  time  delays  determined  from  the 
unsealed  Longshot  traces  was  computed.  Denoted  by  A  in  Figure  6  are  the  results 
for  three  scaled  sets  of  data  with  input  signal-to-noise  ratios  of  lOdB,  5dR  and 
-1.5dB.  The  effects  of  repeating  correlation  after  employing  wavefront  corrections 
are  indicated  by  □  ,  in  Figure  6. 

Typical  waveforms  of  the  scaled  and  unsealed  traces  are  shown  in  Figure  7. 

For  the  latter,  the  signal-to-noise  ratio  is  almost  35dB  and  the  time  delay  can  be 
determined  with  little  difficulty.  However,  when  the  signal-to-noise  ratio  is  reduced 
by  36dB  so  that  the  input  signal-to-noise  ratio  is  negative,  and  even  when  the 
traces  are  appropriately  time-phased  as  is  the  case  in  Figure  7,  it  is  obvious 
that  a  manual  time  delay  measurement  is  impractical. 


II. 1.4  Conclusions 


The  importance  of  accurate  time  delays  cannot  be  over  emphasized.  First,  a 
loss  in  beam  gain  is  experienced^  anc]  the  resulting  beam  defocusing  will  degrade 
array/processing  performance.  Second,  event  location  errors  (see  Appendix  II. 2) 
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b.  Beam  Correlation  Coefficient 


c.  Center  Seismometer  F-4  Subarray 


d.  Beam  Correlation  Coefficient 


Figure  5.  Waveform  Correlation  Traces 
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Figure  6.  Correlation  Estimation  of  Phase  Delays 
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Figure  7  Seismorre'er  Record  cf  Longshot 


are  strongly  dependent  on  the  time  delays.  Third,  the  maintenance  of  both  the  time 
delay  and  location  conversion  libraries  in  a  form  suitable  for  updating  implies  an 
automatic  processing  procedure  for  both  time  delay  determination  and  use. 

The  method  discussed  above  provides  the  desired  results.  The  sparse  experi¬ 
mental  results  available  to  date  warrant  consideration  of  the  technique  as  an  integral 
part  of  event  processing. 


II. 2  EVENT  LOCATION  ACCURACY 


II. 2 . 1  Introduction 


The  error  generated  when  estimating  the  location  of  a  seismic  event  received 
on  a  large  array  can  be  attributed,  in  part,  to  the  presence  of  seismic  noise  at  the 
time  of  arrival  of  the  event  signal  (i.e.,  to  the  fact  that  the  signal-to-noise  ratio 
is  finite),  and,  in  part,  to  imperfect  "calibration."  In  the  present  context  the  term 
"calibration"  refers  to  the  method  and  data  used  to  convert  wavefront  measurements 


(such  as  the  reciprocal  wave  velocities  u 


,  u  )  to  a  geographic  location  for  the  event. 

x  y 


In  Section  2,  an  idealized  model  is  used  to  examine  the  effect  of  seismic  noise 
on  the  accuracy  of  the  estimated  event  location.  In  Section  3,  an  attempt  is  made  to 
incorporate  the  errors  resulting  from  imperfect  calibration. 


II. 2. 2  The  Effect  of  Seismic  Noise  on  Estimated  Event  Location 

Because  of  the  presence  of  seismic  noise  at  the  time  of  arrival  of  a  seismic 
event  at  the  array  site,  errors  are  generated  in  the  estimated  values  of  the  recipro¬ 
cal  phase  speeds  (ux,  Uy)  of  the  event  wavefront.  Hence,  even  with  perfect  "calibra¬ 
tion",  errors  in  the  geographic  event  location  will  result.  These  errors  tend  to 
increase  with  a  decreasing  signal-to-noise  ratio. 

To  assess  the  effect  of  seismic  noise  in  a  quantitative  manner,  we  consider  an 
idealized  large  array  model,  wherein  array  beams  are  formed  from  subarfay  beams 
for  which  the  following  behavior  is  postulated: 

a.  Seismic  noise  received  at  any  one  subarray  is  statistically  independent  from 
that  received  at  any  other  sub  array 

b.  The  signal  waveforms  received  at  all  subarrays  are  identical  except  for 
(frequency  independent)  scale  factors  and  for  travel  time  differences 

c.  The  wavefront  shape  is  accurately  known  and  the  time  delay  for  the  k-th 
subarray  can  be  represented  as  t^  +  xkux  f  Ykuy»  where  t\K  is  the  (accurately 
known)  station  correction  for  the  k-th  subarray  and  (x^,  y^)  are  the  posi¬ 
tion  coordinates  of  the  center  of  that  subarray. 
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The  first  assumption  is  in  reasonable  accord  with  present  experience.  The 
second  assumption  is  not  always  fulfilled:  substantial  decorrelation  has  been 
observed  for  signal  waveforms  received  on  different  subarrays.  However,  by  and 
large,  the  correlation  is  satisfactory  for  a  few  seconds  following  the  signal  onset. 
The  third  assumption  postulates  perfect  "calibration"  of  the  wavefront  shape. 

This  assumption  is  realistic  if  the  seismic  event  originated  in  close  geographic 
proximity  of  previously  recorded  strong  seismic  events.  Where  such  previously 
received  events  are  lacking,  the  assumption  is  not  always  valid. 

With  the  postulates,  the  estimated  reciprocal  phase  speeds  (ux,  uv)  approximate 
Gaussian  variables  with  the  true  quantities  (ux,  Uy)  as  mean  values.  Furthermore, 
assuming  application  of  (frequency  independent)  weighting  of  the  subarray  beams  to 
achieve  maximum  array  beam  output  signal-to-noise  ratio,  the  contours  of  constant 
probability  density  for  (ux,  u,r)  coincide  with  the  constant-loss  contours  of  the  beam 
pattern  centered  at  (ux,  Uy) .  ^Such  weighting  is  desirable,  for  example,  when  the  sub¬ 
arrays  do  not  all  have  the7 same  number  of  seismometers. 

Denoting  by  p  the  probability  for  (ux,  Uy)  to  fall  outside  the  beam  contour  cor¬ 
responding  to  a  loss  of  L  dB  for  a  center  frequency  of  f  Hz,  the  following  approxima¬ 
tion  holds: 


(1) 


where  ctq  is  the  standard  deviation  of  the  timing  error  for  the  signal  waveform  in 
the  LASA  beam  output  and  is  given  by: 


(2) 


Here  S/N  is  the  output  signal-to-noise  ratio  of  the  array  beam  steered  to  the  event 
location  (ux,  Uy),  while  WT  is  the  time -bandwidth  product  (for  LASA,  WT  ~  1  ). 

The  above  formulae  provide  minimum  error  bounds  and  tend  to  under  estimate  the 
errors  in  (ux,  Uy)  particularly  for  low  S/N. 

The  preceding  results  conform  to  the  intuitively  obvious  fact  that  the  location 
error  (ux,  Gy)  can  be  minimized  by  positioning  the  subarrays  in  such  a  manner 
that  the  resulting  beam  pattern  has  minimal  main-lobe  width.  For  the  particular 
case  where  the  subarrays  are  constrained  to  be  located  in  a  circular  area  of 
radius  R,  they  should  be  positioned  on  the  circumference.  In  that  case  formula 
(1)  reduces  to: 


p^exp  -(7rfrR)2  •  ^  •  2WT 


(3) 


where  r  is  the  distance  between  the  estimated  location  (ux,  Uy)  and  the  true  location 
(ux,  Uy).  The  root-mean-square  of  r  is  given  by: 


(4) 
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The  error,  r,  is  measured  in  see/Km  and  is  related  to  the  (ux,  Uy)  plane.  To 
translate  it  as  a  geographic  location  error,  it  must  be  multiplied  by  a  sensitivity 
factor  C  which  generally  depends  on  the  range  from  the  event  location  to  array  as 
well  as  on  the  azimuth  from  the  true  event  location  to  estimated  event  location. 

After  averaging  over  this  azimuth,  the  factor  C  varies  from  about  1.2  x  10^  to  2.0  x 
105  Km^/see.  Thus,  C  =  1 .6  x  10^  Km^/see  appears  to  be  a  satisfactory  compro¬ 
mise  particularly  suitable  for  ranges  from  the  event  location  to  the  array  in  exeess 
of  65°. 

As  an  example,  consider  the  case  where  f  =  1  Hz,  WT  =  1.  Then,  the  root- 
mean-square  geographic  error  is  given  by: 


„„„  ,  .  0.36  x  10  N  'T, 

RMS  geographic  error  «  - — -  J  —  Km,  (5) 

with  the  array  radius  R  being  measured  in  Km.  This  result  is  depicted  graphically 
in  Figure  8. 

We  emphasize  that  minimization  of  the  main- lobe  width  may  not  be  a  praetieal 
design  criterion  since  it  can  lead  to  a  poor  side-lobe  structure. 

The  expressions  (3)  and  (4)  presuppose  that  all  subarrays  show  the  same  output 
signal-to-noise  ratio.  When  this  is  not  the  case,  the  value  of  S/N  in  (3)  and  (4)  must 
be  derated.  For  example,  when  all  subarrays  have  the  same  number  of  seismometers, 
except  for  one  subarray  whieh  has  a  larger  number,  the  value  of  S/N  must  be  derated 
by  less  than  0.4  dB,  provided  that  the  large  subarray  has  no  more  than  half  of  all  the 
seismometers  in  the  array  and  provided  that  the  subarrays  are  suitably  positioned 
relative  to  eaeh  other. 


II. 2. 3  The  Effeet  of  Imperfect  nCalibrationM  on  Location  Error 

The  considerations  in  Section  2  have  assumed  perfect  ’Calibration"  (i.e.,  have 
disregarded  the  errors  arising  when  converting  the  estimated  reciprocal  phase 
speeds  (ux,  fty)  to  geographic  event  location).  These  errors  exist  even  when  the 
event  depth  is  known  or  postulated. 

The  main  eause  for  "calibration"  errors  lies  in  the  fact  that  strong  seismic 
events,  needed  to  establish  the  data  base  for  the  ’Calibration"  technique,  are  rather 
rare;  favor  specific  geographic  areas  to  the  virtual  exclusion  of  other  large  regions; 
have  finite  signal-to-noise  ratio;  and  are  sometimes  difficult  to  pinpoint  in  depth 
and  geographic  location.  These  eircumstanees  not  only  impede  precise  determina¬ 
tion  of  the  relative  station  corrections  required  for  beamforming,  but  also  hamper 
the  establishment  of  an  accurate  conversion  method  from  (u  r,  u  )  to  geographic 
location.  X  V 


Although  extensive  data  analyses  have  been  undertaken  to  obtain  the  relative 
station  corrections  for  LASA,  no  explicit  study  has  yet  been  made  to  calibrate  a 


conversion  method  from  (u 

'  V 


u  )  to  geographic  location. 


Presently,  ray  tracing 
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Figure  8.  Location  Error  with  Perfect  "Calibration" 
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techniques,  based  on  the  Jeffreys-Bullen  tables  and  others,  are  used  for  this 
conversion,  but  no  systematic  effort  has  been  undertaken  to  compare  the  calculated 
event  locations  with  those  obtained  from  other  observations  (such  as  from  a  net¬ 
work  of  seismic  stations)  and  to  establish  corrections  for  possible  bias  in  the  event 
locations  calculated  from  LASA  data.  In  particular,  little  information  appears  to 
exist  regarding  the  event  location  accuracy  presently  achieved  by  LASA. 

One  may  conjecture  that  the  nealibrationM  error  will  eonsist  of  two  portions: 
one  part  being  proportional  to  the  reciprocal  of  the  array  diameter,  2R,  and 
refleeting  the  effect  of  the  finite  array  aperture;  the  second  part  being  independent 
of  aperture  size  and  reflecting  the  location  inaccuracies  for  seismic  events  used  in 
the  data  base.  We  shall  assume  that  both  portions  ean  be  looked  upon  as  random 
errors  when  considered  on  a  worldwide  basis.  Hence,  combining  these  errors  in 
RMS  fashion  with  the  error  presented  by  formula  (4),  we  obtain  the  following  expres¬ 
sion  for  the  total  geographic  RMS  location  error: 


A 


+  B  Km, 


(6) 


+ 


R 


where  the  constants  A  and  B  will  depend  on  the  '’calibration”  accuracy. 

Although  in  the  initial  phases  of  the  operation  of  an  array,  large  "calibration" 
errors  may  exist  (reflected  in  large  values  of  A  and  B  in  formula  (6)),  as  time  pro¬ 
gresses  and  the  database  broadens,  the  "calibration"  errors  should  decrease  and 
the  values  of  A  and  B  should  become  smaller.  Also,  A  and  Bwill  tend  to  be  less 
for  geographic  areas  where  seismic  events  are  relatively  prevalent  than  for  areas 
where  they  are  rare. 

The  results  of  formula  (6)  have  been  plotted  in  Figure  9  using  the  same  data 
as  in  Section  2  (i.e.,  taking  f=l  Hz,  WT=1,  and  Ol.GxlO^  Km^/sec)  and  choosing 
A=10*  Km^  and  B=0.  This  choice  for  A  and  B  implies  a  calibration  error  of  100 
Km  for  an  array  radius  of  100  Km  and  has  been  somewhat  arbitrary. 

It  is  interesting  to  compare  this  error  with  the  location  error  which  would  be 
incurred  when  the  sole  contributing  eause  is  assumed  to  be  the  fact  that  the  sampling 
rate  is  finite,  thus  limiting  the  accuracies  with  which  the  time  delays  for  array 
beamforming  ean  be  implemented.  This  location  error  ean  be  calculated  from 
formula  (1)  provided  that  the  standard  deviation  be  defined  by: 

2  _  (AT)2 

a0  12.  K 

(based  on  the  assumption  that  timing  errors  due  to  time  sampling  are  uniformly 
distributed  over  the  sampling  interval  AT  and  are  uneorrelated  between  subarrays). 
Here,  K  is  the  number  of  subarrays.  For  an  array  whose  subarrays  are  of  equal 
size  and  are  positioned  on  a  circle  with  radius  R,  the  RMS  error  is  found  to  be 
2(70  C/R.  For  a  sampling  rate  of  ten  samples  per  second  (AT  0.1  sec),  K  21  and 
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Figure  9.  Location  Error  with  Imperfect  "Calibration" 


27 


5  2 

C  2x10  Km  see,  one  finds  a  location  error  of  20  Km  when  It  100  Km.  Thus,  we 
assumed  for  Figure  9  a  "calibration"  error  five  times  larger  than  the  error  caused 
by  the  sampling  rate.  For  LASA,  formula  1  leads  to  an  It  MS  error  of  40  Km, 
caused  by  sampling  rate. 


II.3  FIDELITY  OPTIMUM  PROCESSING 


II.3.1  Introduction  and  Summary 

This  appendix  contains  the  results  of  a  theoretical  study  to  estimate  the  ability 
of  fidelity  optimum  processing  to  discriminate  against  one  of  two  simultaneous 
events.  It  compares  the  signal-to-noise  ratio  of  conventional  processing  and  fidelity 
optimum  processing  when  the  noise  consists,  in  part,  of  a  second,  interfering  seis¬ 
mic  event.  Antenna  patterns  of  the  resulting  optimum  processor  are  also  obtained. 

Specifically,  the  noise  is  assumed  to  be  the  sum  of  an  uncorrelated  background 
and  an  interfering  event.  The  background  is  assumed  to  have  the  same  power  spec¬ 
trum  at  each  seismometer.  The  interfering  event  is  also  assumed  to  be  a  second 
order  stochastic  process  with  a  power  spectrum  identical  to  that  of  the  background 
noise. 

In  this  situation,  the  gain  from  optimum  processing  is  a  function  of  the  array 
configuration,  the  radius  of  the  array,  the  ratio  of  the  power  of  the  interfering 
event  to  the  background  noise  power,  and  the  magnitude  of  the  difference  in  inverse 
phase  velocity  between  the  event  of  interest  and  the  interfering  event.  Three  array 
configurations  were  considered:  The  LASA  configuration,  21  seismometers  equally 
spaced  on  a  circle  and  five  seismometers  equally  spaced  on  a  circle.  If  we  let  R 
be  the  radius  of  the  array  and  |Au  |  be  the  magnitude  of  the  difference  in  inverse  phase 
velocity  between  the  two  events,  then  the  gain  is  a  continuous  function  of  R  |  Au|  and 
we  determined  this  function  over  the  range  of  interest.  Ratios  of  interfering  event 
noise  to  background  noise  varying  from  0.05  to  100  were  considered. 

Substantial  gains  were  found  to  occur  for  r|au  |  >  0.2  for  the  LASA  configuration 
and  for  R|  Au|  >0.15  for  the  other  configurations.  In  the  above  ranges  of  R  |au|  and 
for  large  ratios  of  interfering  event  noise  power  to  background  noise  power,  the 
optimum  processor  essentially  eliminated  the  interfering  event. 


II.3.2  The  Power  Output  of  a  Linear  Processor 

We  suppose  we  have  K  seismometers  and  designate  their  sampled  outputs  by 
X.(t),...,  X^(t).  We  assume: 

Xk(t)  =  S(t)  +  Nk(t),  k  =  1 . K;  (8) 

where  each  seismometer  output  is  the  sum  of  a  signal,  which  we  are  interested  in 
estimating,  and  noise.  The  signals  are  assumed  to  be  lined  up  or  in  phase  from 
seismometer  to  seismometer.  We  suppose  the  noises,  N^(t),  to  form  a  multi- 
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variate  second  order  process  with  a  spectral  matrix  P^,(f).  That  is, 


E{Nk(t)  Nk,(t+r)}  =  <^k,(r) 


and 

*  (9) 

Pkk'  (f)  At  £,  \k'  (T)  e  ^  ’ 

7*=_oo 

where  At  =  sampling  interval. 


Now  we  will  consider  the  case  where  a  linear  operator  with  impulse  response 
0k(s)  and  transfer  function  Ak  (f)  is  applied  to  the  kth  seismometer  output  for  k=l, 

. .  .,K.  These  filter  outputs  are  summed  to  obtain  the  output  of  the  array,  Y(t).  This 
is  illustrated  in  Figure  10. 


0k(s)  and  Ak(f)  are  related  by: 


=  ^  At(t)e2l)fs  df, 


Vs> 


-w 


and 


Afc(f)  =  At  0k(s)  e 


-  2  7T j  1 S 


S=-00 


where  w  - 


2  At 


(10) 


We  will  be  interested  in  the  output  noise  power  of  an  array  with  such  processing. 
The  output  noise  power  spectrum  (see  Figure  11)  is  given  by: 

1  V)  Ak,(f)  pkk,(f).  (11) 

k^; 

By  Ak(f),  we  mean  the  complex  conjugate  of  Ak(f). 

Hence,  the  output  noise  power  is  given  by: 
w 

I  I  Ak'f>Ak'«t>pkk'(f> df-  <12> 

"w  k,k' 

II.3.3  Form  of  the  Fidelity  Optimum  Processor 

In  fidelity  optimum  processing,  the  goal  is  to  minimize  the  output  noise  power 

subject  to  the  constraint  that  the  signal  is  undistorted.  Now,  the  transfer  function 

K 

seen  by  the  signal  is  Z  A,  (f);  hence,  the  goal  of  fidelity  optimum  processing  is  to 

k=l  K 

minimize: 
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Figure  10.  Linear  Processor 
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S  1  V 


co  Ak-(f>  pkk.  <'> 

-w  k,k' 

subject  to  the  condition: 


(i:t) 


(f)  =  i. 


(14) 


7  2  8 

It  can  be  shown  (see  Capon  or  Vanderkulk  ’  )  that  the  Aj_(f)  satisfying  the  above 
are  given  by: 


Ipkk’ 


<f) 


Ak(f)  - 


£  Pki> 

k,k' 


(15) 


where  P.  ,  ,  (f)  is  the  k,k'  element  of  the  inverse  of  the  KxK  matrix  P,  ,  .  (f) . 
k,k'  '  kk' '  ' 

II.3.4  Inverting  the  Spectral  Matrix 

Pkki  (f)  for  fixed  f  is  a  matrix  of  complex  numbers.  If  we  let: 

Pkk-  <f>  -  Ckk'  ([>  +  J  «V  <*>  (16) 

and 

pkk'(f>  =  Dkk'<f> +  J  Ekk'<f>- 

then  it  is  easily  seen  that  the  following  relationship  holds  between  real  2K  x  2K  matrices 


-Q 


D  -E\ 


Q 


E 


D 


(17) 


/ 


,-l 


This  relationship  can  be  used. to  find  P,  ,  ,  (f)  with  a  real  matrix  inversion  routine. 
We  can  also  use  the  above  form  to  show  that 

Pkk'  (f)  =  pDo.  (*)  and’  hence’  that: 


k'k 
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.-1 


kk'  |  kk'  (  l_,  kk' 

k,k'  k,k'  ^  J  k,k' 

C  -Q  \ 

Consider  first  the  matrix  |  I  .  Since  C  is  symmetric  and  Q  is  skew 

Q  Cj 


(18) 


D 


symmetric,  this  matrix  is  symmetric.  Hence,  its  inverse 


,  is  sym- 


D 


metric  which  means  that  E  =  -E. 

By  E7  we  mean  the  transpose  of  E.  This  is  equivalent  to  P  (f)  =  P^q.  (f)  • 


II. 3. 5  A  Convenient  Form  for  the  Output  Noise  Power  with  Optimum  Processing 
We  saw  earlier  that  the  output  noise  power  in  optimum  processing  was  given 


by: 


w 


^  ‘  ^  Ak(f)  Ak'(f)  Pkk'(f)  df' 


(19) 


-w  kk' 

This  is  an  integral  over  a  quadratic  form  which  can  be  written  as: 

I  V>  Ak'  <f>  Pkk-  <f>  =  A  T  P  A. 
kk' 


(20) 


and  P  is  the  cross  spectral  matrix.  We  have  suppressed  f  for  simplicity.  However, 
we  also  saw  that: 


>  pki> 

T' 


AJf)  = 


pkk-  ® 


=  P_1  A  , 


(21) 


kk' 


where 


A  = 


Pkk'  (f) 


(\\ 

i  / 


kk' 
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T  -1  -1 

/"«w»  *  X  /-«W»  T  _  I  rvj  T  TT 

Hence,  A  PA=AP  P  P  A  =  A  P  A 

-1  ^~1~^  ~  ~T  ^  -1 

However,  since  P^,  (f)  =  Pttt  (f),  we  have  A  =  A  and  P  =  P 


kk' 


Hence,  finally  : 


~T  T  -1  kk1 

A  P  A  =  A  P  A  =  — 


I  Pkk'(f| 


Pkk'(t> 


kk' 


(T  pkk'<f>''_1 

\  kk1 


(22) 


or 


A  P  A  = 


This  yields  • 


-1  \ 
Re{pkk-  (f» 


kk' 


-1 


/ 


W 

output  noise  power  =  II  Ak(f)  Ak.(f)  pkk.(f) 


-w  kk. 


(23) 


(24) 


or 


w 


-II  Re{pkk-i(t>i 


-w 


kk' 


df. 


(25) 


II.3.6  The  Output  Noise  Power  in  Conventional  Processing 

In  conventional  processing  the  output  noise  power  is  given  by: 

U/ 

(1/K2)  ^  Re{Pkk,(f)}  df 


I 


-w 


kk' 


II. 3. 7  The  Special  Case  of  Noise  Consisting  of  a  Single  Seismic  Disturbance  Plus 
an  Uncorrelated  Background 

Suppose  the  noise  nk(t)  to  be  the  sum  of  two  noises  n^  k(t)  and  n  k(t).  Further 
suppose  nc  k(t)  to  be  a  seismic  disturbance  with  delays  t  .  ,r  relative  to  the 
delays  of  S(t).  Then:  K 


ncfkW  = 


Rc^-Tk 


). 


(20) 
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Further  suppose  that  nc  (t)  is  a  second  order  process  with  power  spectrum  (f) 
and  autocovariance  function  (p  (T)-  Then: 

E{nck(t)  nck(t  +  T)l  =  E(nc (t_Tk)  nc(t_Tkt  +  T)^  (27) 


=  <P(t+  Tk  -  rk,) 

Taking  the  Fourier  transform  we  obtain  Pc  (f)  e  J  '  k  k!/  as  the  contribution 
to  the  cross  spectrum  P^k,  (f) . 

Now  njj  k(t)  is  assumed  to  be  a  background  noise  which  is  uncorrelated  from 
seismometer  to  seismometer  and  uncorrelated  with  the  seismic  disturbance.  If 
it  has  a  power  spectrum  P^  (f)  then  the  cross  spectra  are  given  by 

Pkk.(t)  =  Pb(f)  <5  (k>k’)  +  pc(f)  e2?rjf(Tk"  Tk’);  (28) 

or  if  we  let  P  (f)  /  P  (f)  =  M(f) } 

0  u 

Pkk'(f)  =  pb  (f)  {  6  (k-k’)  +  M(f)  e2?rjf(Tk  "  V},  (29) 

8  2  ‘j 

where  <5  (kkf)  is  the  Kronecker  delta.  From  Vanderkulk  and  Kelly  and  Levin  *  we  have 


(f) 


pb(f) 


{6  (kk’)  -  M  (f)  e27rjf(Tk  "  Tk!)  /  (1  +  KM  (f))}  . 


(30) 


Using  the  above  we  can  write  down  expressions  for  the  transfer  functions  A^(f) 
and  the  noise  power  output  for  optimal  and  conventional  processing.  The 
transfer  function  is: 


Ak(f) 


Ml 

l+KM  (f) 

M  (f) 

l+KM  (f) . 


e  -2?rjfTk' 


27Tjf  (\ 


Tw) 


kk' 


(31) 


The  output  noise 


(optimum  processing)  is: 


M(l) 

l+KM(f) 


I 


e  27T  j f  ( Tk  -  V 


df, 


(32) 
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and  the  output  noise  power  (conventional  processing)  is: 


w 


f  Pb  (f)  {K+M  (f)  ^  e27r^f(Tk  "  V}  /  K2  df. 


(33) 


kk' 


The  calculations  can  be  simplified  by  noticing  that: 


2 


e2')f<Tk  "V 


27t jf  x  1 2 

e  kl 


(34) 


kk' 
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II. 3. 8  Determination  of  Delays 

We  suppose  K  seismometers  are  located  by  the  vectors  p^  with  polar  coordinates 
(Rrk,  T^) .  R  is  a  parameter  which  controls  the  overall  radius  of  the  array.  We  let 
u”l  be  the  inverse  phase  velocity^ofJ,he  event  of  interest,"^  be  the  inverse  phase  veloci¬ 
ty  of  the  interfering  event  and  Au=u1-U2*  Further  we  let  Au  have  polar  coordinates 
(  |au  |,  <p) .  With  these  definitions  the  delays,  t,  ,  of  the  interfering  event  with  respect 
to  the  event  of  interest  are  given  by 


Tk  =  Au  .  pk 


(35) 


=  R  |Au  |  rk  cos  (yk  -  <p) 

In  our  investigation  the  will  be  normalized  so  that  their  maximum  is  approximately 
equal  to  one.  Hence,  the  delays  are  a  function  of  R  |Au(  and  <p. 

II.3.9  The  Antenna  Pattern 

For  the  class  of  linearly  processed  arrays  under  consideration,  the  antenna 
gain  in  the  direction  with  delays  r^,  . .  . ,  is  : 


I  Ak  <f> 


2 


antenna  gain  (f,  t^,  .  .  .  ,  r. ) 


(36) 


k 
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This  is  a  power  gain,  and  if  we  have  an  input  from  this  direction  with  power 
spectrum  P(f),  the  power  output  will  be: 


w 


(37) 


I  p,f>  I 


-w 


k 


II.3.10  Specific  Calculations 

We  now  give  the  results  of  some  calculations  to  determine  the  ability  of  optimum 
processing  to  discriminate  against  the  noise  described  in  Section  II. 3. 7  (i.e.,  noise  con¬ 
sisting  of  a  single  seismic  disturbance  plus  an  uncorrelated  background). 

In  these  calculation,  we  assumed  Pc  (f)  and  P|j(f)  to  have  the  same  shape:  i.e., 

Pc(f)  =  P(f)=MPb(f)-  We  took  as  the  shape  of  P(f)  a  smooth  approximation  to  the 
spectrum  of  the  output  of  an  individual  seismometer  for  the  Kamchatka  earthquake 
(April  8,  1966)  as  given  in  Figure  3-43,  Page  3-130  of  Ref.  4. 

As  the  ratio  M  of  the  spectra  Pc  (f)  and  P|3(f),we  took  values  from  0.05  to  100.  We 
considered  three  configurations  of  seismometers.  The  first  was  that  for  the  21 
subarrays  of  the  present  I.ASA.  The  r^  andy^  are  given  in  Table  1.  The  second 
was  an  array  of  21  seismometers  located  equally  spaced  on  the  circumference  of 
a  circle.  The  third  was  an  array  of  five  seismometers  located  equally  spaced  on  the 
circumference  of  a  circle. 

For  these  configurations  and  for  a  wide  range  of  directions  of  arrival  of  the 
unwanted  seismic  event,  we  obtained  the  noise  power  output  for  optimum  processing 
and  for  conventional  processing.  We  also  computed  the  ratio  of  the  conventional 
noise  power  output  to  the  optimum  noise  power  output.  In  addition,  we  calculated 
antenna  patterns  for  both  a  sample  case  of  optimum  processing  and  conventional 
processing. 

In  these  calculations,  we  assumed  that  the  array  was  pointed  in  the  direction 
corresponding  to  all  delays  equal  to  zero.  The  noise  event  will  then  have  a  set  of 
delays,  not  all  equal  to  zero,  which  we  saw  in  Section  II. 3. 8  are  functions  of  It  |Au|  and 
0.  where  R  is  the  radius  of  the  array  and  (Au,  $)  are  polar  coordinates  of  the  difference 
in  inverse  phase  velocity  of  the  event  of  interest  and  the  interfering  event.  We  did 
the  calculations  for  various  values  of  <t>  and  generally  for  a  range  of  It  |Au|  >_  0.  The 
results  are  plotted  against  R  |au  |  with  <p fixed. 

II. 3. 10.1  Results  of  the  I.ASA  Configuration 

In  Figure  12,  we  have  plotted  the  noise  power  output  for  conventional  processing 
and,  in  Figure  13,  the  noise  power  output  for  optimum  processing  as  a  function  of 
the  direction  of  the  interfering  seismic  event.  This  is  done  for  values  of  M  from 
0.05  to  100.  Specifically,  the  noise  spectral  matrix  is  assumed  to  be: 
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Table  1 


LAS  A  COORDINATES 


Instrument  Group 

Normalized  Position 

rk 

AO 

0.000 

B1 

0.123 

55.2° 

B2 

0.007 

141.7° 

B3 

0.082 

246.1° 

B4 

0.090 

346.7° 

Cl 

0.183 

23.8° 

C2 

0.163 

97.8° 

C3 

0.131 

192.5° 

C4 

0.128 

293.7° 

D1 

0.305 

56.8° 

D2 

0.263 

141.9° 

D3 

0.252 

231.9° 

D4 

0.308 

336.1° 

El 

0.542 

13.6° 

E2 

0.686 

106.9° 

E3 

0.607 

183.3° 

E4 

0.538 

278.3° 

FI 

1.093 

46.4° 

F2 

1.037 

147.4° 

F3 

1.036 

219.4° 

F4 

0.973 

235.4° 
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Figure  12.  Output  Noise  Power,  Conventional  Processing, 
LASA 
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Figure  13.  Output  Noise  Power,  Optimum  Processing,  LASA 
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Pkk,(f)  =  P(f)  {6  (kk')  +  Me 


2rif  (rk  -  rk,)} 


(38) 


where  the  rk  and  Tk,  are  determined  by  the  direction  of  the  interfering  seismic 
event  relative  to  the  event  of  interest  as  described  in  Section  II. 3. 8.  We  have  assumed 
0  =  0  (i.e.,  all  events  are  arriving  from  the  same  relative  azimuthal  direction)  and 
varied  the  r|Au|  for  different  values  of  M. 

In  Figure  14,  we  have  plotted  the  dB  gain  in  noise  rejection  of  optimum  process¬ 
ing  over  conventional  processing  for  the  LASA  configuration  and  values  of  M  from 
1  to  100.  Notice  that  substantial  gains  occur  for  r|au|>0.2. 

In  Figure  15,  we  have  plotted  the  antenna  gain  of  the  optimum  processor  for  an 
event  with  R|au|  =  0.3,  0=0  and  M  =  10.  Cuts  radially  out  from  the  center  are 
plotted  for  =  0,  n/8,  . . . ,  n  .  Notice  the  sharp  null  in  the  direction  of  the  interfering 
event.  Figure  16  gives  several  radial  cuts  of  the  conventional  antenna  pattern. 

II. 3. 10. 2  Results  for  21  Seismometers  Located  Equally  spaced  on  a  Circle 

Figure  17  contains  the  noise  power  output  for  conventional  processing  and 
Figure  18  the  noise  output  power  for  optimal  processing  for  the  case  for  21  seismo¬ 
meters  equally  spaced  on  a  circle  for  M  from  0.05  to  10.  In  Figure  19,  we  have 
plotted  the  dB  gain  in  noise  rejection  from  optimum  processing  for  M  from  1  to  10. 
There  is  a  significant  improvement  over  the  21  seismometer  configuration  of  the 
LASA  array.  Figure  20  contains  the  radial  cut  of  the  optimum  processor  antenna 
pattern  through  the  interfering  event  for  the  same  interfering  event  depicted  in 
Figure  15. 


II. 3. 10. 3  Results  for  Five  Seismometers  Located  Equally  Spaced  on  a  Circle 

Figures  21,  22,  23,  and  24  give  some  results  for  a  configuration  of  five  seismo¬ 
meters  located  equally  spaced  on  a  circle.  Figure  23  gives  the  antenna  pattern  of 
the  optimum  processor  again  for  the  same  interfering  event  as  in  Figures  15  and  20. 
Figure  24  gives  the  antenna  pattern  of  the  conventional  processor  for  0 <0<7r  /5.  The 
remainder  of  this  pattern  can  be  obtained  by  symmetry. 

II. 3. 11  Optimum  Processing  Consisting  of  a  Single  Weighting  of  Each  Seismometer 


Output  Before  Addition 


Suppose  we  wish  to  restrict  our  processing  to  simple  weighted  sum  of  the  seis¬ 
mometer  outputs.  That  is,  we  wish  our  output  to  be  of  the  form: 


(.-!<)) 


k 
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(9P)  woQ 
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Figure  14,  dB  Gain  Due  to  Optimum  Processing,  LASA  Configurati 


Antenna  Pattern  (dB) 


15 


Interfering  R|au  |-  0.3 

Event  <p  =  0 

M  =  10 


Figure  15.  Antenna  Pattern 

Optimum  Processing 
Interfering  Event 
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Antenna  Pattern  (dB) 


Figure  16.  Antenna  Pattern 

Conventional  Processing 
LASA  Configuration 
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Figure  17.  Output  Noise  Power 

Conventional  Processing 
21  Seismometers  on  a  Circle 
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Output  Noise  Power  (dB) 


20 


Figure  18.  Output  Noise  Power 
Optimum  Processing 
21  Seismometers  on  a  Circle 
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(gp)  uiDQ 
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Figure  19.  dB  Gain  Due  to  Optimum  Processing 
21  Seismometers  on  a  Circle 


Antenna  Pattern  (dB) 


r|au| 


Figure  20.  Antenna  Pattern,  Optimum  Processing 
21  Seismometers  on  a  Circle 
Interfering  Event 
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Figure  21 .  Output  Noise  Power 

5  Seismometers  on  a  Circle 


49 


(BP)  u!dO 
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Antenna  Pattern  (dB) 


10 


-10 


-20 


-30 


0.15  0.30  0.45  0.60  0.75  0.90 

R|au| 

Figure  23.  Antenna  Pattern 

Optimum  Processing 
5  Seismometers  on  a  Circle 
Interfering  Event 
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10 


Figure  24.  Antenna  Gain 

Conventional  Processing 
5  Seismometers  on  a  Circle 
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with  the  fidelity  constraint: 


l 


(40) 


This  is  the  same  as  restricting  the  (f)  of  sections  II. 3. 2  and  II. 3.3  to  be  of  the  form 
Ak(f)  =  cfc.  With  this  restriction  the  output  noise  power  is  given  by: 

w 


output  noise  power 


I  I  ck  V  pkk'  <*> df 

-w  k,k' 


(41) 


w  w 

^  Ck  Ck'  Pkk'  (f)  df 


kk' 


-w 


w 


-  y  °k  ck'  y  Re{pkk'  <f>) df- 


kk' 


-w 


(42) 


(43) 


Now  if  we  let: 
w 

$  Re<pkk-(t)l  df  =  “kk-  <44) 

-w 


Then  we  wish  to  minimize 


I 


kk' 


Ck  Ck' 


II 

kk' 


(45) 


subject  to  the  constraint: 


I 


(46) 
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This  is  the  same  abstract  problem  described  in  Section  II. 3. 3  and  the  optimum  set  of 
c^  are  given  by: 


(47) 


Further  as  in  Section  It. 3. 5: 
output  noise  power 


■  I 


Ck  V  11 


k,k' 


kk' 


-1 


(48) 


(49) 


and  from  Section  II. 3. 9  the  antenna  gain  is: 
antenna  gain  (opt.  processing)  = 

k 


ck  e-2lJfTk|2. 


(50) 


For  the  case  of  present  interest,  a  noise  consisting  of  an  interfering  event 
with  spectrum  P  (f)  and  delays  r,  and  a  uncorrelated  background  with  spectrum 
(f) ,  we  have  C 

v  - 1  pkk''f» 

=  6  (k,k')  J  Pb(f)  df  +  j*  Pc(f)  cos  27Tf(Tk  -  TRt)  df  ; 
or  if  Pb(f)  =  P(f)  =  Pc(f)  /M, 


(51) 


nkk,  =  6(k,k')  J  P(f)  df  +  M  J  P(f)  cos  27rf(rk  -  r^) 


df. 


(52) 


We  made  some  calculations  comparing  the  gain  in  noise  rejection  using  this 
simple  processing  with  the  gain  from  complete  optimum  processing.  Figures  25 
and  26  give  the  noise  power  output  for  the  LASA  configuration  and  for  21  seismo¬ 
meters  on  a  circle  for  single  multiplication  optimum  processing.  There  is  con¬ 
siderable  improvement  over  conventional  processing  but  substantially  less  than 
with  complete  optimum  processing.  For  five  seismometers  on  a  circle,  a  sin¬ 
gle  calculation  for  the  value  M=5  showed  single  multiplication  optimum  process¬ 
ing  to  give  almost  no  improvement  over  conventional  processing. 
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Output  Noise  Power  (dB) 


30 


Figure  25 


Output  Noise  Power,  Single  Multiplication 
Optimum  Processing  LASA  Configuration 


55 


25 


0.15  0.30  0.45  0.60  0.75 

Figure  26.  Output  Noise  Power 

Single  Multiplication  Optimum  Processing 
2  1  Seismometers  on  a  Circle 


0. 

R| 


II. 4 .  FIXED  POINT  FFT  EXECUTION 


II.4.1  Introduction 


II. 4. 1.1  Summary 

This  appendix  contains  the  results  of  an  effort  directed  at  estimating  the  impact 
of  microcoding  on  the  performance  of  the  fast  Fourier  transform  algorithm  (FFT) 
on  the  IBM  System/360  Model  40.  There  were  two  major  aspects  to  this  effort; 
estimating  the  speed  of  a  fixed  point.  Model  40,  microcode  and  estimating  the  error 
in  the  fixed  point  calculation.  The  microcode  and  its  speed  are  discussed  in  detail 
in  Section  II. 4. 3. 

If  N  is  the  number  of  complex  points  to  be  transformed  (N  is  assumed  to  be  a 
power  of  two),  then  the  following  estimate  is  obtained  for  the  speed  of  a  fixed  point 
microcode  operating  on  16  bit  numbers; 

Transform  Time  ~  52  N  logi;>N  ji  secs. 

This  is  approximately  ten  times  faster  than  both  floating  and  fixed  point  machine 
language  codes  for  the  algorithm  on  the  Model  40. 

The  accuracy  of  the  fixed  point  calculation  is  discussed,  in  detail,  in  Section  EL 4. 2. 
The  fixed  point  calculation  we  considered,  both  for  microcoding  and  in  the  accuracy 
estimation,  involved  keeping  the  array  properly  scaled  by  rescaling  at  any  stage  in 
which  there  occurred  an  overflow.  With  this  scheme,  the  error  is  a  function  of  the 
numbers  and  timing  of  these  overflows.  However,  an  upper  bound  to  the  error  can 
be  obtained  by  assuming  an  overflow  and  rescaling  at  each  stage  of  the  calculation. 
Again,  if  we  assume  we  are  transforming  N  =  2^  complex  points,  then  the  upper 
bound  of  the  ratio  of  the  root  mean  square  error  to  the  root  mean  square  of  the 
resulting  transform  increases  as  the  square  root  of  N  or  by  vlT  at  each  of  the  M 
stages  of  the  calculation.  Specifically  for  the  Model  40  microcode  which  operates 
on  16  bit  numbers,  we  have: 


RMS  (error) 
RMS  (result)  - 


2(M+3)/2  -15  |U  3, 

RMS  (initial  array) 


This  is  well  within  the  accuracy  requirements  of  the  FFT  applications  envisaged  for 
LASA  data.  This  point  is  discussed,  in  detail,  in  subsection  EI.2.4.5. 


II.4.1.2  The  Finite  Fourier  Transform 

If  X(j)  j  =  0,  1,  .  .  .  ,  N-l  is  a  sequence  of  complex  numbers  then  the  finite 
Fourier  transform  of  X(j)  is  the  sequence: 
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N-l 

A(m)  =  (1/N)  V  X(j)e"27rljm/N 


m=0,  1,  .  .  .  ,  N-l. 


(53) 


J=° 


The  inverse  transform  is 

X(i)  =  £  A,m)e21I^m//N 

m=0 


(54) 


In  both  of  the  above  equations  i  =  (-1)  .  We  will  be  considering  a  fixed  point 

calculation  of  these  transforms  using  the  Fast  Fourier  Transform  algorithm9  ,  10  . 
In  connection  with  Equation  (53) ,  we  will  actually  obtain  N  A  (m)  from  X  (j).  We  will 
then  give  this  result  with  N“1  as  part  of  an  overall  scale  factor.  Considering  the 
calculation  of  N  A  (m)  from  X  (j)  or  X  (j)  from  A  (m)  Parseval's  theorem  states: 


N-l  N-l 

y  N)i 2  =  n  > 

1  A(m)  |2 

j-o  m=o 

N-l 

N-l 

or  \  |  NA(m)  |2  =  N 

I  lx«> 

L — i 

m=o 

j~° 

(55) 


and  we  see  that  the  mean  square  value  of  the  result  is  N  times  the  mean  square 
value  of  the  initial  sequence.  This  fact  will  be  used  below. 


II. 4. 1.3  The  Inner  Loop  of  the  Fast  Fourier  Transform  Algorithm;  Step  by  Step 
Scaling 

The  inner  loop  of  the  FFT  algorithm  operates  on  two  complex  numbers  from 
the  array.  It  takes  these  two  numbers  and  produces  two  new  complex  numbers 
which  replace  the  original  ones  in  the  array.  Let  Xm(i)  and  Xm(j)  be  the  original 
complex  numbers.  Then,  the  new  pairXm+  ^(i),  ^m  +  ^(j)  are  given  by: 


X  ^  (i)  =  X  (i)  +  X  (j)W 
m+1'  ’  m  mu/ 

X  ^  (j)  =  X  (i)  -  X  (j)\V 
m+1  'J/  m'  ’  m 
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where  W  is  a  complex  root  of  unity.  If  we  write  these  equations  out  in  terms  of 
their  real  and  imaginary  parts,  we  get: 

Re{Xm+i(i)}  =  Re{Xm(i)}  +  Re{Xm(i)}  Re{W}  -  Im{Xm(i)}  Im{w} 

(57) 

Im^Xm+l(j)^  =  Im(Xm(i)}  "  Re^Xm(^  IrTW  "  Im^-Xm(j^  ReM 

At  each  stage  the  algorithm  goe^ihrough  the  entire  array  of  N  numbers  in  this 
fashion,  two  at  a  time,  If  N  =  2  ,  then  the  number  of  such  stages  in  the  computation 

is  M. 

As  we  move  from  stage  to  stage  through  the  calculation,  the  magnitudes  of  the 
numbers  in  the  array  generally  increase  which  means  that  the  array  can  be  kept 
properly  scaled  by  overflow  tests  and  right  shifts.  Consider  first  the  root  mean 
square  of  the  complex  numbers.  From  equations  (56)  we  have 


Hence,  in  the  root  mean  square  sense,  the  numbers  (both  real  and  complex)  are 
increasing  by  V~2  at  each  stage.  Consider  next  the  maximum  modulus  of  the  com¬ 
plex  numbers.  From  equations  (56)  one  can  easily  show  that 


mux  {|xm(i)|,|xm(j)|}  <  max  {|Xm+1(i)|,  |Xm+1(j)|}  <.  2  max  {|  Xm(i)  |  ,  |xm<j)|} 

(59) 


Hence,  the  maximum  modulus  of  the  array  of  complex  numbers  is  non-decreasing. 

Both  of  the  above  results  suggest  that  the  fixed  point  calculation  be  done  as 
follows;  the  entire  original  array  and  the  Wfs  are  scaled  with  the  numbers  moved  as 
far  to  the  left  as  possible;  after  the  multiplication,  the  most  significant  bits  are 
retained;  after  the  addition,  there  is  a  test  for  overflow;  if  overflow  is  detected,  the 
entire  array  is  shifted  right  one  bit  and  the  calculation  is  continued.  Equation  (59) 
shows  that  the  maximum  modulus  of  the  complex  numbers  cannot  increase  by  more 
than  a  factor  of  2.  Hence,  if  scaling  were  controlled  by  this  modulus,  there  would 
be  no  more  than  one  shift  per  stage.  However,  the  method  we  have  outlined  above 
would  control  the  scaling  on  the  basis  of  overflow  of  the  calculation  of  the  real  and 
imaginary  parts.  It  can  be  shown  that,  in  this  case,  it  is  possible  (although  unlikely) 
to  get  more  than  one  overflow  per  stage  or  a  single  overflow  of  more  than  one  bit. 

It  is  impossible  to  get  an  overflow  of  more  than  two  bits.  (This  can  be  seen  from 
equations  (57)  since  there  are  only  two  additions). 
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11.4.2  An  Accuracy  Analysis  of  the  Fixed  Point  Calculation 


II  .4. 2.1  Introduction 

We  will  assume,  in  this  analysis,  that  the  inputs  (i.e.,  the  real  and  imaginary 
parts  of  X(j)  or  A(m))  are  represented  by  B  bits  plus  a  sign.  We  assume  the  binary 
point  lies  to  the  left  of  the  leftmost  bit.  We  showed  earlier  that  the  magnitudes  of  these 
members  of  the  array  would  generally  increase  as  we  moved  from  stage  to  stage  in  the 
calculation.  Hence,  the  method  of  operation  is  to  test  for  overflow  within  the  inner 
loop.  If  there  is  no  overflow,  the  calculation  proceeds  as  usual.  If  there  is  an  over¬ 
flow,  then  the  two  inputs  producing  the  overflow  are  shifted  right  until  there  is  no 
overflow.  The  amount  of  the  shift  is  recorded  (it  will  be  either  one  or  two  bits) 
and  the  entire  array  is  shifted  right  this  same  amount.  In  this  scheme,  wc  shift 
not  only  those  elements  we've  already  calculated  but  also  those  yet  to  be  done. 

(This  latter  shifting  could  be  done  in  subsequent  calculations  if  this  were  more 
efficient).  The  total  number  of  shifts  is  accumulated  and  the  power  of  two,  raised 
to  the  negative  of  this  total  number  of  shifts,  constitutes  an  overall  scale  factor 
to  be  applied  to  the  final  array. 

Our  numbers  in  most  applications  will  come  from  an  A-D  converter  and  will 
contain  a  quantization  error.  We  will  not  consider  the  propagation  of  this  quantiza¬ 
tion  error.  Our  purpose  is  to  compare  the  accuracy  of  a  short  word  fixed  point 
computation  with  a  longer  word  floating  point  computation.  In  practice,  the  quanti¬ 
zation  error  would  be  common  to  both  alternatives. 

There  are  two  operations  which  produce  errors  which  are  propagated  through 
the  calculation. 

a.  When  two  B  bit  numbers  are  multiplied  together  a  2B  bit  product  results. 

If  this  product  is  rounded  to  B  bits,  an  error  whose  variance  is: 

a2  _  2~2B/12  (60) 

is  created.  This  error  has  a  standard  deviation: 

A  =  2"B//l2  ~0.3(2_B)  (61) 

b.  When  two  B  bit  numbers  are  added  together  and  there  is  an  overflow, 
then  the  sum  must  be  shifted  right  and  a  bit  lost.  If  this  bit  is  a  zero, 
there  is  no  error.  If  it  is  a  one,  there  is  an  error  of  +2"®  depending 
upon  whether  the  number  is  positive  or  negative.  The  variance  of  this 
error  (it  is  unbiased  assuming  there  are  an  equal  number  of  positive  and 
negative  numbers)  is: 

a2  =  2-2B/2>  (62) 
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It  has  a  standard  deviation 


A2  =  2~B  1  « 0.7(2  B) 


(63) 


11.4.2.2  Upper  Bound  Analysis 


In  this  section,  we  give  an  upper  bound  analysis  of  the  ratio  of  the  RMS 
error  to  the  R  M  S  of  the  answer.  This  upper  bound  is  obtained  by  assuming 
that  during  each  step  of  the  calculation  there  is  an  overflow  and  a  need  to  rescale. 
We  let  X|<  be  a  typical  real  element  at  the  kth  stage  (i.e.,  the  real  or  imaginary 
part  of  a  complex  element)  and  let 


V(XjJ  =  variance  of  X^. 


2  2 

We  will,  in  this  analysis,  replace  A^  by  GA 


r 


We  will  also  let 


(64) 


Since  the  first  stage  gives  an  overflow,  the  original  data  must  be  rescaled 
or  truncated  by  one  bit.  Hence, 


V(X0)  =  6A2.  (65) 

In  going  from  the  original  data  to  the  results  of  the  first  stage,  W=1  and,  hence  there 
is  no  multiplication  and  we  either  add  or  subtract.  Further,  we  assume  that  the  next 
stage  will  result  in  an  overflow  and  hence  we  will  have  to  rescale.  This  gives 


=  2V(XQ)  +  4  •  6A2 
V(X])  -  2(6 A2)  +  4  •  6 A2. 


(66) 


In  these  expressions  and  this  entire  discussion,  we  are  assuming  all  errors  to  be 
independent  and,  hence,  that  the  variance  of  the  sum  is  the  sum  of  the  variances. 
Going  from  the  first  stage  to  the  second  stage,  we  have  W  =  (-1)  ■  ^  and  again  there 
are  only  additions  and  subtractions.  Thus,  with  the  rescaling 


V(X.,)  -  2V(X1)  +  42  •  6A2 

-  22(6A2)  +  2(4- 6A2)  +  42  •  6A2. 


(67) 


In  going  from  the  second  stage  to  the  third  stage,  we  have  multiplications  and  wc 
have  them  in  all  subsequent  stages.  In  generating  the  third  stage,  half  the  inner 
loops  have  multiplications.  Consider  the  first  equation  of  equations  (57).  All  the 
other  equations  are  identical  in  terms  of  error  propagation.  Remember  that  X^  (i) 
is  complex; 
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(68) 


Re  |X3(i)}  =  Re  {X2  (i)}  +  Re{X2  (j)}  Re{W}  -  Im{X2(j)}  Im{W} 


Equation  (68)  yields,  with  rounding  to  B  bits  after  the  addition  and  with  rescaling. 


V'(X3)  =  V(X2)  +  [Re2{X2(j)}  +  Im2{X2(j)}]  V(W) 

+  [Re2(W)  +  Im2(W)]  V(X0) 


(69) 


+  (43A2)  +  43  •  6A2 

=  V(X2)  +  |X2(j)|2A2  +  V(X2)  +  (43A2)  +  4 3  •  6A2 


In  Equation  (69),  the  first  term  is  the  variance  of  the  first  term  of  (68).  The  second 
and  third  terms  of  (69)  are  the  variance  of  the  full  2B  bit  products  given  by  the 
second  and  third  terms  of  (68).  The  fourth  term  of  (69)  is  the  result  of  rounding 
after  the  addition.  The  fifth  term  is  the  rescaling  term.  Finally,  we  saw  in  Equa¬ 
tion  (58)  that  the  average  modulus  squared  of  the  complex  numbers  is  increasing 
by  a  factor  of  2  every  stage.  Hence,  if  we  let  K  equal  the  average  modulus  squared 
of  the  initial  array,  i.e., 


N-l 


i=0 

then  we  have: 


V'  (X3)  =  2V(X2)  +  22KA2  +  (4 3 A2)  +  43  •  6  •  A2 


(70) 


Equation  (70)  would  be  correct  for  V(X^)  if  all  the  inner  loops  involved  multiplica¬ 
tions.  However,  at  this  stage  only  half* of  them  do  and  hence: 


In  the  next  stage,  three  quarters  of  the  inner  loops  require  multiplications  and 
these  multiplications  get  progressively  more  numerous  as  the  stages  increase. 
Hence,  from  here  on,  we  will  assume  all  stages  have  multiplications  in  all  the  inner 
loops.  Thus,  applying  the  above  techniques,  we  get: 
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V(X4)  =  24(6A2)  +  23(6  •  4A2)  +  22  (6  •  42A2)  +  2(6  •  43A2)  +  6  •  44A2 
+  22KA2  +  23KA2  +  43A2  +  44A2; 


(72) 


and  generally  if  M  is  the  last  stage: 

V(XM)=  2M(6A2)  +  2M_1(6-4A2)  +  .  .  .+  2(6-4M_1a2) 

+  2M+2KA2  +  (M-3)2M_1KA2  +  2M_4(43A2) 
+  2M_4(44A2)  +  .  .  .+  (4Ma2) 


=  (1.5)2M  +  2A2(l+2+ .  .  .+2M  *)  +  (M-2.5)2M  *KA 


0M+3  2  0M+4  _  0M-4  2 

+  2  A  +  2  (1+2+ .  .  . +2  )  a  • 


(73) 


Or 


V(XM)  «  (1.5)22M+2A2  +  (M-2.5)2M  *KA2  +  2M+3A2  +  22M+1a2 
^22M+3A2  +  (M-2.5)2M-1KA2  +  2M+3A2- 


(74) 


K  is  the  average  of  the  square  of  the  absolute  values  of  the  initial  complex  array. 
Hence,  applying  Parseval’s  theorem,  Equation  (55),  the  average  of  the  square  of 
the  absolute  values  of  the  final  array  will  be  2^K.  What  is  most  meaningful  in 
this  case,  however,  is  the  mean  square  of  the  real  numbers  which  is  2^K/2. 
Hence,  we  have: 


V(X  )  _  2M+3A2  (M-2 .5)  A2/2  23A2 

- —  ~  -  +  -  +  - 

\T  * 

2mK/2  K/2  1/2  K/2 

and  finally, 


RMS  (error) 

RMS  (result 


M+3 


K/2 


M+3 

2  2  2~B(0.3) 

RMS  (initial  array) 


(75) 


(76) 
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Equation  (76)  gives  an  approximate  upper  bound  for  the  ratio  of  the  RMS 
of  the  error  to  R  M.S.  of  the  answer.  Notice  that  this  bound  increases  as  the  \rN 
or  1/2  bit  per  stage. 


II.4.2.3  Lower  Bound  Analysis 

We  will  now  obtain  a  lower  bound  for  the  ratio  of  R  M  S  of  the  error  to  the 
R  M  S  of  the  answer.  We  obtain  this  lower  bound  by  assuming  that  there  are  no 
overflows  in  the  calculation  and,  hence,  no  shifts  of  the  array.  In  this  case, 


V(X0)  =  V(X1)  =  V(X2)  =  0.  (77) 

In  the  third  stage,  half  of  the  inner  loops  involve  a  multiplication  and  hence: 

V(X3)  =  (1/2)  (22KA2)  +  1/2  (A2).  (78) 

This  can  be  seen  by  considering  the  first  term  of  equations  (69).  The  first  term 
of  (78)  comes  from  the  second  term  of  the  first  of  equations  (69).  The  second  term 
of  (78)  is  caused  by  the  rounding  to  B  bits.  Now,  as  before: 


V(X4)  =  2V(X3)  +  23KA2  +  A2 

=  22KA2  +  23KA2  +  A2  +  A2. 


(79) 


Finally, 


V(XM)  .  2M-2KA2  ♦  (M-3)2M-1KA2  .  2M-V  +  2M‘5A2 

_M-b  .2  .2 

+  2  A  +  .  .  .  +  A 

=  (M-2.5)2M_1KA2  +  2M_3A2  +  (1+.  .  ,+2M_5)  A2 


(80) 


~  /*/r  o  Cw,M_1TrA2  0M-3.2  0M-4  .2 

~  (M-2.5)2  KA  +2  A  +  2  A  . 

As  in  subsection  2.2,  the  mean  square  of  the  array  of  resulting  real  numbers  is 
2M  .  k/2.  Hence,  we  have: 


V*XM> 

2MK/2 


(M-2.5)  a 


A^/S 

K/2 


2/c 
A  /6 

K/2 


(81) 


(M-2.5)  A  . 
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and 


RMS  (error) 

RMS  (result) 


(82) 


(M-2.5) 


1/2 


The  lower  bound  increases  as  M  =  log2N.  This  is  the  rate  of  increase  which 
has  been  observed  for  the  floating  point  calculation.  ’ 1 


II.4.2.4  Some  Experimental  Results 

An  IBM  7094  program  was  written  to  perform  a  fixed  point  calculation  using 
the  fast  Fourier  transform  algorithm,  as  described  above.  The  program  was 
capable  of  simulating  a  fixed  point  machine  of  any  word  size  up  to  35  bits  plus  a 
sign.  Experiments  were  run  with  fixed  point  numbers  of  17  bits  plus  a  sign.  This 
corresponds  to  B  =  17  in  the  analysis  of  subsections  II. 4. 2. 2  and  H.4.2.3. 

The  experiments  were  performed  in  the  following  way.  Floating  point  input 
was  fixed  to  17  bits  plus  a  sign.  This  fixed  input  was  then  transformed  with  the 
fixed  point  program.  The  fixed  point  output  was  then  floated.  Next,  the  fixed 
point  17-bit  input  was  floated  and  a  floating  point  transform  taken.  Since  this 
floating  point  transform  uses  a  floating  point  word  with  a  27-bit  mantissa,  it  was 
considered  the  correct  answer.  Finally,  the  R  M  S  of  the  difference  between  the 
fixed  point  and  floating  point  answers  was  taken.  We  also  obtained  the  maximum 
absolute  error  and  average  error. 

Figure  27  contains  the  result  ol  transforming  random  numbers  which  lie  between 
zero  and  one  (placed  in  both  the  real  and  imaginary  parts).  In  this  and  subsequent 
tests,  three  runs  were  made  for  every  power  of  2  from  8  to  2048.  Since  these  ran¬ 
dom  numbers  have  a  DC  component  of  one-half,  the  fixed  point  program  must 
rescale  at  least  every  stage  but  one.  Hence,  one  would  expect  the  error  to  lie  close 
to  the  theoretical  upper  bound  as  given  by  Equation  (76).  This  theoretical  upper 
bound  is  also  plotted  in  Figure  27  and  the  results  are  seen  to  lie  slightly  above  it. 

The  R  M  S  of  the  original  array,  K/2,  is  approximately  0.6. 

Figure  28  contains  the  results  of  transforming  three  sine  waves  plus  random 
numbers  between  zero  and  one-half  in  the  real  part  and  all  zeros  in  the  imaginary 
part.  Specifically: 

Re{X(j)}  -  1/2 [Y (j)  +  (1/2)  sin  (27r8j/N) 

+  (1/4)  sin  (27r4j/N)  +  (1/4)  sin  (27T8j/N) J 

Im{X(j)}  -  0 
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Figure  27.  Experimental  Error  Results:  Random  Numbers 
Between  0  and  1,  B=17 


(4|fis3y)swa/(JOJJ3)SWy 
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Figure  28.  Experimental  Error  Results:  Random  Numbers  Plus  3  Sine  Waves 
0  <  Random  Numbers  <  1/2,  B  =  17 


where  the  Y(j)  are  random  numbers  between  zero  and  one.  Again,  there  is  a 
DC  component  of  magnitude  one-fourth  and  the  array  must  be  rescaled  at  least 
at  every  stage  but  two.  Thus,  one  would  expect  these  results  to  be  lower  relative 
to  the  theoretical  upper  bound  than  the  case  depicted  in  Figure  27.  From  Figure  28 
one  can  see  that  this  is  in  fact  the  case.  The  R  M  S  of  the  original  array,  K/2, 
is,  in  this  case,  approximately  0.35.  This  is  the  reason  the  upper  bound  curve  is 
higher  than  that  of  Figure  27. 

Figure  29  contains  the  results  of  transforming  random  numbers  from  minus 
one  to  one  (in  both  real  and  imaginary  parts).  In  this  case,  the  DC  component  is 
zero  and  there  is  no  other  strong  component.  The  number  of  shifts  should  be 
approximately  (log2N)/2  or  one-half  shift  per  stage.  Hence,  one  would  expect 
the  error  curve  to  lie  well  below  the  theoretical  upper  bound  as  is  the  case.  In 
this  case,  K/2  =  0.6. 

Figure  30  contains  the  results  of  an  experiment  identical  to  that  corresponding 
to  Figure  28,  except  the  random  numbers  are  between  plus  and  minus  one  half.  The 
results  are  as  expected.  In  this  case,  K/2  =  0.35. 

Finally,  Figure  31  contains  the  results  of  transforming  a  sine  wave  in  the  real 
part  and  zero  in  the  imaginary  part.  The  sine  wave  was  sin  (27rj/8).  Although  in 
this  case  the  array  must  be  rescaled  in  at  least  every  stage  but  two,  the  error  is 
well  below  the  upper  bound.  Here,  K/2  =  0.5. 

These  calculations  were  made  with  rounding  rather  than  truncation.  The  error 
with  truncation  was  somewhat  larger.  Rounding  can  be  included  in  the  microcode 
described  in  Section  II.4.3  with  a  negligible  decrease  in  speed.  The  bias,  as 
reflected  by  the  average  error,  was  in  general  negligible  compared  with  the 
RMS  error.  The  maximum  error  was  roughly  what  would  be  expected. 


II.4.2.5  Conclusions  and  Additional  Comments 

Figure  32  below  gives  the  upper  bound  on  the  ratio  of  R  M  S  of  the  error  to  the 
R  M  S  of  the  answer  assuming  B  =  15  and  assumming  that  the  RMS  of  the  initial 
array  is  three-tenths  (i.e.,  K/2  =  0.3).  This  corresponds  to  the  microcoding 

situation  for  the  Model  40. 

There  are  two  main  uses  to  which  the  fast  Fourier  transform  would  be  put  in 
seismic  signal  processing.  The  first  would  be  in  estimating  the  spectral  matrix 
of  a  subarray  of  seismometer  outputs  as  a  first  step  in  the  calculation  of  the  coef¬ 
ficients  for  optimum  processing  (see  13,  14,  2  or  8).  The  second  would  be  in 
estimating  a  sequence  of  short  term  spectra  of  the  LASA  beam  output  as  one  of 
the  first  steps  in  discrimination  processing.  In  both  of  these  cases,  an  acceptable 
and  computationally  efficient  procedure  is  to  take  transforms  of  short,  overlapping 
segments  of  the  time  series;  calculate  periodograms  or  cross-periodograms  and 
time  average  over  these  to  increase  the  statistical  accuracy.  (See  Section  8  of 
ref.  11,  and  ref.  13.) 
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Figure  29.  Experimental  Error  Results:  Random  Numbers 
Between  -1  and  1,  B=17 
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Figure  30.  Experimental  Error  Results:  Random  Numbers  Plus  3  Sine  Waves 
-  1/2  <  Random  Numbers  <  1/2,  B  =  17 


CO  so 

I 

o 


CM 


o 


NO 


CM 


LO 

I 

o 


(4|nssy)sWH/(JOJJ3)SWil 


l\ 

II 

CO 


71 


Figure  31.  Experimental  Error  Results:  Single  Sine  Wave, 
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Figure  32,  Error  Upper  Bound  for  15-Bit  Calculation 


As  concerns  the  fixed  point  accuracy,  these  remarks  are  important  in  two 
respects.  First,  since  short  segments  are  transformed,  the  fixed  point  error 
would  be  kept  small.  In  the  spectral  matrix  estimation,  the  size  of  the  segments 
would  be  of  the  order  of  the  number  of  coefficients  in  the  optimum  processing 
procedure.  In  the  estimation  of  LASA  beam  spectra,  it  would  be  of  the  order  of 
the  reciprocal  of  the  resolution  desired  (in  terms  of  the  Nyquist  frequency  as 
unity).  In  both  of  these  cases,  256  would  be  a  high  figure  for  the  length  of  the 
segments  to  be  transformed.  For  segments  of  length  256,  the  fixed  point  error 
would  be  small  compared  to  the  statistical  error  in  the  spectral  estimates.  Its 
only  effect  might  be  a  slight  one  on  the  available  dynamic  range.  Second,  the 
time  averaging  would  reduce  this  fixed  point  calculation  error  as  well  as  the 
statistical  error.  Hence,  for  the  purposes  of  spectral  estimation  important  in 
LASA  signal  processing,  a  fixed  point  calculation  appears  to  be  more  than  ade¬ 
quate. 

As  a  final  comment,  since  the  magnitude  of  the  error  has  a  sensitive  and  direct 
relationship  to  the  number  of  times  the  array  is  scaled,  it  would  be  best  to  remove 
any  substantial  DC  component  which  happened  to  be  present. 
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II. 4. 3  Speed  Advantage  of  the  Microcoded  FFT 


II. 4. 3.1  The  Inner  Loops 

M 

The  FFT  algorithm  operates  on  an  array  of  2  complex  numbers,  taking  M 
stages  on  the  array.  Each  stage  comprises  2^-1  passes  through  the  inner  loop 
computations,  where  each  pass  processes  two  complex  points  of  the  array.  In 
terms  of  these  points,  called  OPerand  X  and  OPerand  Y,  and  the  complex  root 
of  unity  U,  the  inner  loop  is  either  a  Type  I  or  Type  II: 


Type  1: 


Complex 

OPX 

OPY 

=  OPX  +  OPY 
=  OPX  -  OPY 

Real 

ROPX 

ROPY 

=  ROPX  +  ROPY 
=  ROPX  -  ROPY 

Imaginary 

IOPX 

IOPY 

=  IOPX  +  IOPY 
=  IOPX  -  IOPY 

Type  II: 

Complex 

OPX 

OPY 

=  OPX  +  OPY •  U 
=  OPX  -  OPY • U 

Real 

ROPX 

ROPY 

=  ROPX  +  (ROPY  -  RU  -  IOPY-  IU) 

=  ROPX  -  (ROPY  -  RU  -  IOPY-  IU) 

Imaginary 

IOPX 

IOPY 

=  IOPX  +  (ROPY  -  IU  +  IOPY-  RU) 

=  IOPX  -  (ROPY-  IU  +  IOPY-  RU) 

The  Type  I  inner  loop  is  used  for  all  passes  in  the  first  and  second  stages 
through  the  array.  Thereafter,  its  application  is  halved  for  each  subsequent 
stage  (see  Table  2). 

The  number  of  Type  II  inner  loop  executions  in  the  total  FFT  process  does 
not  exceed  the  Type  I  inner  loops  unless  M  >  5.  Actually,  the  ratio  Type  11/ 
Type  I  is  asymptotic  to  M/2,  but  this  is  not  a  good  measure  of  the  relative  fre¬ 
quencies  of  the  types  unless  M  is  large.  Nevertheless,  the  Type  II  inner  loop 
is  the  dominant  factor  in  the  FFT  process,  because  it  involves  four  real  multi¬ 
plications  and  six  real  additions.  Type  I  involves  no  multiplications,  and  only 
four  additions.  For  this  reason,  our  estimate  has  centered  around  the  Type  II 
inner  loop. 
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Table  2 
INNER  LOOP 


Stage 

Number  of  Type  I 

Inner  Loops 

Number  of  Type  II 

Inner  Loops 

1 

2M-1 

0 

2 

2M-1 

0 

3 

2M-2 

2M-2 

M-l 

22 

02  3  .  „M-2 

M 

2 

2  +22  +  .  .  .  +  2M_2 

II. 4. 3. 2  Type  II  Inner  Loop  Process  and  Flow  Chart 

M 

To  maximize  and  preserve  significance  in  the  fixed  point  operands,  the  2" 
complex  numbers  are  normalized  in  the  sense  that: 

a.  The  high  order  bit  (bit  position  0)  is  a  sign  bit 

b.  The  binary  point  is  assumed  to  be  between  bit  positions  0  and  1 

c.  At  least  one  real  or  imaginary  component  of  one  of  the  complex 
points  has  the  value  1  in  bit  position  1. 

With  this  representation  of  the  complex  points,  the  execution  of  inner  loop 
passes  may  give  rise  to  an  overflow  of  one  or  two  bit  positions.  Actually,  a  two 
position  overflow  has  very  low  probability;  and  a  three  position  or  more  overflow 
is  impossible. 

M-l 

Consider  the  2  inner  loop  passes  in  any  stage.  If  an  overflow  of  one  bit 
position  occurs  first  in  pass  K,  this  overflow  divides  the  stage  into  two  parts: 

a.  Previously  executed  passes;  that  is,  passes  1  to  K-l 

M-l 

b.  Subsequently  executed  passes;  that  is,  passes  K+l  to  2 

Similarly,  previous  and  subsequent  passes  may  be  defined  with  respect  to  a 
two  position  overflow.  This  will  divide  the  stage  into  three  parts,  if  the  two- 
overflow  initially  occurs  in  a  pass  subsequent  to  the  one -overflow. 

When  an  overflow  of  one  or  two  positions  occurs,  the  entire  array  must  be 
shifted  right  one  or  two  bit  positions  to  preserve  significance.  This  is  most 
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expeditiously  done  by  shifting  the  results  of  each  subsequently  executed  pass  in 
the  stage  before  these  results  are  returned  to  main  storage;  and,  after  the  stage 
is  completed,  by  returning  to  the  beginning  and  shifting  the  results  of  all  pre¬ 
viously  executed  passes.  The  latter  process,  of  shifting  previously  executed 
results,  will  be  called  post  shifting.  An  array  is  post  shifted  until  an  array  ele¬ 
ment  with  a  stop-address  is  detected.  This  is  the  address  of  an  inner  loop 
operand  in  the  pass  which  first  gave  rise  to  the  overflow.  Saving  these  addresses 
is  a  bookkeeping  operation  which  must  be  performed  when  an  initial  overflow  of 
one  or  two  places  occurs. 

From  the  above,  we  see  that  the  nature  of  the  inner  loop  process  depends  on 
whether  an  overflow  condition  has  occurred  in  a  previously  executed  pass,  and 
on  whether  a  new  or  additional  shift  condition  has  occurred  in  the  currently  exe¬ 
cuted  inner  loop  pass.  All  of  the  possibilities  are  given  in  Table  3. 

For  each  case  in  Table  3,  we  give  a  probability.  These  assigned  probabili¬ 
ties  are  very  rough  estimates.  They  are  not  to  be  taken  literally,  except  as  an 
indication  that,  for  purposes  of  estimating  microprogram  execution  time,  it 
suffices  to  consider  only  Cases  A,  D,  and  E.  The  effect  of  the  low  probability 
Cases  B,  C,  F,  and  G  is  negligible.  These  conclusions  are  consistent  with  error 
analysis  results. 


Table  3 

INNER  LOOP  OVERFLOW  CONDITIONS 


Case 

Ovfl  Positions 
Previous  Pass 

Ovfl  Positions 
Current  Pass 

Probability 

Total  Shift 
Before  Storing 

A 

0 

0 

0.500 

0 

B 

0 

1 

0.010 

1 

C 

0 

2 

0.001 

2 

D 

1 

0 

0.230 

1 

E 

1 

1 

0.230 

1 

F 

1 

2 

0.001 

2 

G 

2 

0 

0.028 

2 

If  two  16-bit  normalized  binary  numbers  (binary  points  assumed  between  bit 
positions  0  and  1)  are  multiplied  together,  as  with  SUMP  hardware,  the  high  half¬ 
word  of  the  result  will  have  its  binary  point  between  bit  positions  1  and  2.  Hence, 
the  high  halfword  of  the  product  cannot  be  added  to  a  normalized  number  without 
first  either  shifting  the  product  left  one  position,  or  shifting  the  normalized  number 
right  one  position,  to  line  up  the  binary  points.  The  alternative  chosen  depends  on 
the  part  of  the  inner  loop  being  processed  (computation  of  RX  and  RY,  or  computation 
of  IX  and  IY),  and  on  whether  or  not  there  is  a  pending  shift  from  a  previous  pass. 
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A  flowchart  of  the  inner  loop  process,  showing  the  testing  and  handling  of 
the  various  overflow  conditions,  is  shown  in  figures  33,  34,  and  35.  The  states 
Y2  and  Y3  are  used  to  signal  overflow  conditions  from  previously  executed 
passes  as  illustrated  in  Table  4. 


Table  4 

OVERFLOW  CONDITIONS 


Y2 

Y3 

0 

0 

No  Pending  Overflow 

0 

1 

Pending  Overflow  of  1  Position 

1 

1 

Pending  Overflow  of  2  Positions 

In  the  currently  executed  pass,  the  occurrence  of  overflows  in  the  computa¬ 
tion  of  RX  and  IX  is  recorded  in  the  B1  REG,  and  overflows  in  the  computation 
of  IlY  and  IY  are  recorded  in  the  BO  REG. 

New  overflows,  conditions  B,  C,  and  F,  can  be  recognized  either  in  the  cal¬ 
culation  of  RX  and  RY,  or  IX  and  IY.  Thus,  we  have  two  exits  each  for  these  con¬ 
ditions  in  the  flow  charts. 

A  three -character  notation  near  the  lower  right  of  the  flow  chart  boxes  cor¬ 
relates  the  flow  chart  with  the  microprogram  given  on  figures  36,  37,  38,  and  39. 
The  first  character  gives  the  sheet  number,  and  the  second  two  give  the  block 
coordinates  of  the  microinstruction  on  the  sheet. 


II. 4. 3. 3  Execution  Time  Estimates 

The  microcoded  execution  time  estimates  for  the  Type  II  inner  loop  cases 
shown  in  figures  34  and  35  are  given  in  Table  5  and  assume  both  fixed  point 
operands  and  the  SUMP  hardware  for  multiplication. 

Figures  40  and  41  give  the  OS/360  program  for  the  inner  loops.  The  execu¬ 
tion  for  these  is  as  follows: 

Fixed  Point,  assuming  no  pending  shift:  485  psccs 

Fixed  Point,  assuming  one  pending  shift:  480  psecs 

Floating  Point:  560  psccs 
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Figure  33.  Flow  Chart  for  inner  Loop  of  Fast  Fouiier  T-an  form-f- 
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9 

Concurrent  Processing - *■ 


SUMP(RY  •  IU  +  IY  •  RU) 

STORE  RX  and  RY  in  MS 

2J5 

2L1 

1 

(RY  •  IU  +  IY 

•  RU)  LEFT  SHIFT 

3C5 


CASE  A  EXIT  3N3 


Figure  34.  Flow  Chart  for  Inner  Loop  of  Fast  Fourier 
Transform— II 
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4C2 

4C4 

4C8 


CASE  B 


CASES  D,  E 


Figure  35.  Flow  Chart  for  Inner  Loop  of  Fast  Fourier  Transform— III 
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READ  RY  ROPY  NE&ATIVE 
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Figure  36a  inner  Loop  Microprogram  ! 
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F-gure  37„  Inner  Loop  M'croprogram  il 
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F:gure  38,  nner  Loop  M  c  opiog'am  HI 
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F>gure  39,  Inner  Loop  Microprogram  IV 


Table  5 


TYPE  II  INNER  LOOP  EXECUTION  ESTIMATES 


Case 

Probability 

Microinstructions 

Microseconds 

A 

0.5 

86 

54 

B 

0.01 

78 

+  (5)* 

52 

C 

0.001 

81 

+  (20)* 

63 

D 

0.23 

78 

49 

E 

0.23 

78 

49 

F 

0.001 

73 

+  (15)* 

55 

G 

0.028 

(85)* 

53 

Average 

52 

*In  the  cases  not  microcoded,  the  numbers  in  parentheses  are  estimates 
of  the  instructions  to  accomplish  the  following: 

B:  Set  Y3,  save  array  stop-addresses  for  post  shift  processing 

C:  Save  stop  addresses,  shift  operands  2  places,  store 

F:  Save  stop  addresses,  shift  operands  1  place,  store 

G:  Inner -loop  process,  shifting  all  results  two  places. 
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REGISTER  1  CONTAINS  ADDRESS  OF  OPERAND  X,  AOPX 
REGISTER  2  CONTAINS  ADDRESS  OF  OPERAND  Y,  AOPY 
REGISTER  3  CONTAINS  REAL  PART  OF  W,RW,  LEFT  ADJUSTED 
IN  LOW  HALF  WORD 

REGISTER  4  CONTAINS  IW,  LEFT  ADJUSTED  IN  LOW  HALF  WORD 
PROGRAM  1  FIXED  POINT  DATA 


LR 

5,3 

REG5=RW 

LR 

6,4 

REG6=IW 

LR 

7,3 

LR 

8,4 

MH 

5,0(2) 

REG5=RW::RY 

MH 

6,2(2) 

REG6=IW”IY 

MH 

8,0(2) 

REG8=IW”RY 

TM 

PR I OR, X1 

01'  IF  BIT  7  OF  PRIOR  IS  1  CONDITION  CODE  IS  1 

BC 

1,  LET 

BRANCH  TO  LET  IF  AT  LEAST  1  PENDING  SHIFT 

NO  PRIOR  PENDING  SHIFT 

SR 

5,6 

RW::RY-IW::I  Y 

SLA 

5,1 

BC 

l,OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

AR 

7,8 

RW::I  Y+IW::RY 

SLA 

7,1 

BC 

1 , OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

SRA 

5,16 

SRA 

7,16 

LH 

9,0(1) 

REG9-OPERAND  X  REAL  PART 

LR 

10,9 

REGIO^OPERAND  X 

SR 

9,5 

RY=RX-(RW::RY-I  W::I  Y) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

AR 

10,5 

RX=RX+(RW::RY- 1 W « I Y  ) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

STH 

9,0(2) 

STORE  RY 

STH 

10,0(1) 

STORE  RX 

LH 

9,2(1) 

REG 9= I MAG I NARY  PART  OF  OPAND  X 

LR 

10,9 

SR 

9,7 

I  Y=IX-(RW"IY+IW::RY) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

AR 

10,7 

I  X=  I X+  (RW::  I Y+ 1  W::RY  ) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT  PROGRAMMED 

STH 

9,2(2) 

STORE  IY 

STH 

10,2(1) 

STORE  IX 

END  OF  NO  PRIOR  PENDING  SHIFT  CASE 


Figure  40.  OS/360  Programs  and  Execution  Estimates  for 
FFT  Inner  Loop  (I,  II,  and  III) 
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AT  LEAST  ONE  PRIOR  SHIFT  PENDING 


LET 


SR 

5,6 

RW::RY-IW::I  Y 

AR 

7,8 

RW::I  Y+IW::RY 

TM 

PRIOR, X' 

02'  IF  BIT  6  OF  PRIOR 

IS  1  CONDITION  CODE 

BC 

1, ANTHER 

2  PRIOR  SHIFTS, CASE  NOT  PROGRAMMED 

ONE  PRIOR  SHIFT 

SRA 

5,16 

SRA 

7,16 

LH 

9,0(1) 

REG9=RX 

SRA 

9,1 

LR 

10,9 

REG10=RX 

SR 

9,5 

RY=RX-(RW::RY-IW::IY) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT 

PROGRAMMED 

AR 

10,5 

RX=RX+(RW::RY- 1  W::  I Y  ) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT 

PROGRAMMED 

STH 

9,0(2) 

STORE  RY 

STH 

10,0(1) 

STORE  RX 

LH 

9,2(1) 

REG9=IX 

SRA 

9,1 

LR 

10,9 

REG 10= IX 

SR 

9,7 

IY=IX-(RW::IY+IW::RY) 

BC 

l,OVFL 

OVERFLOW  CASE  NOT 

PROGRAMMED 

AR 

10,7 

IX=IX+(RW"IY+IW::RY  ) 

BC 

1 , OVFL 

OVERFLOW  CASE  NOT 

PROGRAMMED 

STH 

9,2(2) 

STORE  IY 

STH 

10,2(1) 

.STORE  IX 

END  OF  ONE  PRIOR  SHIFT  CASE 


Figure  40.  (continued) 
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TIME  FOR  INNER  LOOP  FIXED  POINT 
ASSUMING  NO  PENDING  SHIFT 
INSTRUCTION  NUMBER  TIME  TOTAL  TIMES 


LR 

6 

7.3 

45.0 

MH 

4 

43.0 

180.0 

TM 

2 

9.38 

18.76 

BC 

6 

9.38 

56.28 

SR 

3 

7.5 

22.5 

AR 

3 

7.5 

22.5 

SRA 

16 

2 

18.75 

37.50 

SLA 

1 

4 

10. 

40. 

LH 

2 

18.13 

36.26 

STH 

2 

10.63 

21.26 

TOTALS 

34 

480.06  MICROSECONDS 

TIME  FOR  INNER  LOOP  FIXED  POINT 
ASSUMING  ONE  PENDING  SHIFT 


INSTRUCTION 

NUMBER 

TIME 

TOTAL  TIMES 

LR 

6 

7.5 

45.0 

MH 

4 

45.0 

180.0 

TM 

2 

9.38 

18.76 

BC 

6 

9.38 

56.28 

SR 

3 

7.5 

22.5 

AR 

3 

7.5 

22.5 

SRA  16 

2 

18.75 

37.50 

STH 

4 

10. 

40. 

SRA  1 

2 

18.13 

36.26 

LH 

2 

10.63 

21.26 

TOTALS 

34 

480.06  MICROSECONDS 

Figure  40.  (concluded) 
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REG1  CONTAINS  AOPX, IMAGINARY  PART  AT  AOPX+4 
REG2  CONTAINS  AOPY, IMAGINARY  PART  AT  AOPY+4 
FLREGO  CONTAINS  RW 
FLREG2  CONTAINS  IW 


LER 

4,0 

FLREG  4  =  RW 

LER 

6,2 

FLREG  6  =  IW 

ME 

4,0(2) 

RW”RY 

ME 

6,4(2) 

IW”IY 

SER 

4,6 

RW"RY-IW”IY 

LE 

6,0(1) 

RX 

SER 

6,4 

RX-(RW::RY-IW”IY) 

STE 

6,  TEMP 

TY  TEMP  STORE 

LE 

6,0(1) 

RX 

AER 

6,4 

RX+(RW”RY-IW”IY) 

STE 

6,0(1) 

RX  STORED 

LER 

4,0 

RW 

LER 

6,2 

IW 

ME 

4,4(2) 

RW::I  Y 

ME 

6,0(2) 

IW-RY 

AER 

4,6 

RW”IY  +  IW”RY 

LE 

6,4(1) 

IX 

SER 

6,4 

IX-(RW”IY+IW”RY) 

STE 

4,4(2) 

STORE  IY 

LE 

6,4(1) 

IX 

AER 

6,4 

IX  +  (RW”IY+IW”RY) 

STE 

6,4(1) 

STORE  IX 

MVC 

4(4,1), 

TEMP 

STORE  RY 

END  OF 

FLOATING  POINT  PROGRAM 

TIME  FOR 

INNER  LOOP 

FLOATING  POINT 

INSTRUCTION  NUMBER 

TIME 

TOTAL  TIMES 

LFR 

4 

7.5 

30.0 

ME 

4 

80.63  322.52 

SER 

3 

14.3 

42.9 

LE 

4 

11.88  45.76 

STE 

4 

12.5 

50.0 

AER 

3 

14.3 

42.9 

MVC 

1 

26.25  26.25 

TOTALS 

23 

560.33  MICROSECONDS 

Figure  41.  Floating  Point  Data 
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These  figures,  divided  by  the  52  /isecs  microcode  estimate,  give: 


Advantage  over  OS/360  Fixed  Point:  9.2 

Advantage  over  OS/360  Floating  Point:  10.8 


Remarks  on  Above  Results 

a.  Note  that  the  inner  loop  with  0-pending  and  0-found  overflow  is  longer 
in  S/360  execution  than  the  case  with  1 -pending  and  0-found.  This 

is  so  because  the  normalized  format  of  the  data  gives  a  natural  right 
shift  of  one  place  in  multiplication.  This  must  be  unshifted  in  the 
no-pending  shift  case;  in  the  1-pending  shift  case  the  result  is  cor¬ 
rectly  positioned  after  the  multiply,  without  further  processing. 

Taking  advantage  of  this  gives  rise  to  the  discrepancy  between  Case  A, 
and  cases  D  and  E,  making  the  pending  shift  cases  D  and  E  faster. 

b.  The  floating  point  S/360  process  does  not  appear  unduly  longer  than 
the  fixed  point  (560  to  480).  Note  that  the  floating  point  process 
involves  no  testing  for  overflow,  and  no  special  provisions  for  shifting. 

c.  In  the  microcoded  cases  a  pending  shift  is  processed  by  shifting  the 
result  of  the  innerloop  computation.  Alternatively,  a  pending  shift  may 
be  executed  on  the  input  data  before  the  inner  loop  processing. 

If  input  data  is  shifted,  there  is  an  improvement  in  speed  over  the  output 
result  shifting  because  some  ALU  idle  time  during  SUMP  operation 
can  be  utilized  for  the  process.  But  the  preshift  loses  accuracy: 

Note  that  if  ct  denotes  shift, 


ctx  +  ay  =  <r(x+y) 

For  o  ( 01)  +  a  (01)  =  0  +  0  =  0 
but  ct  (01  +  01)  =  o-(10)  =  1 

The  amount  of  time  saving  in  executing  the  pending  shift  prior  to  inner 
loop  execution  will  be  about  five  microinstructions,  or  about  five  per¬ 
cent  of  the  inner  loop. 

d.  The  estimates  given  are  for  the  power-of-two  FFT  process. 
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11.4.3.4  The  Overall  FFT  Process 


Besides  the  inner  loop,  the  Fast  Fourier  Transform  comprises: 

a.  Index  Processing,  for  accessing  operands  for  the  inner  loop 

b.  Initialization,  instruction  fetch  and  setting  initial  counter  values 
and  working  local  storage. 

c.  Post  shifting,  suppose  an  operand  array  has  256  data  points.  The 
inner  loop  is  executed  128  times  in  each  stage.  If  an  overflow  occurs 
first  in  the  64th  execution,  all  inner  loop  outputs  from  the  64th  to  the 
128  will  be  properly  positioned  in  storage.  But  the  stage  will  not  be 
completed  until  the  outputs  from  the  first  63  inner  loop  executions 
are  also  right  shifted,  in  the  post  shifting  process. 

d.  Data  scramble  routine,  for  bit  reversal  of  operand  addresses.  A 
FORTRAN  program  (Figure  42),  flowchart  (Figure  43)  and  microcode 
estimate  (Figure  44)  of  the  scramble  routine  are  given.  The  micro¬ 
programming  is  illustrated  by  figures  45  and  46.  It  is  interesting  to 
note  that  the  advantage  for  the  scramble  routine  is  appreciably  higher 
than  the  usual  seven  or  eight  for  a  scientific  application.  In  general, 
the  advantage  of  microcode  derives  from  the  elimination  of  I-Fetches 
in  a  sequence  of  system  instructions,  and  from  the  determinate  oper¬ 
ands.  In  the  microcoded  scramble  routine,  there  is  another  factor: 

The  "branch  on  condition"  is  virtually  free  in  microcode,  but  costs  a 
system  instruction.  The  scramble  routine's  high  advantage  results 
from  the  highly  branching  nature  of  the  routine. 

It  is  well  known  that  the  number  of  inner  loop  executions  in  the  total  FFT 
processing  of  2^  =  N  complex  points  is  M  •  2^-l  (M  is  the  number  of  stages 
through  the  array,  and  2  points  of  an  array  are  processed  with  each  inner  loop 
execution).  Experience  indicates  that  the  overall  FFT  execution  time  may  be 
estimated  as  twice  the  inner  loops  execution  time: 

2  •  M  •  2M_1  -  t  •  N  •  log0N, 

where  the  initial  2  provides  for  items  a  through  d  above,  and  where  t  is  execu¬ 
tion  time  per  inner  loop.  Our  estimate  of  t  is  52  microseconds,  so  that: 


52 N  log^N,  usees 

is  the  microcoded  FFT  process  time  for  an  array  of  N  points.  A  few  values  are 
given  in  Table  6. 
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Figure  42.  Jennrich  Scramble  Fortran  Program 
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Figure  43.  Jennrich  Scramble  Flow  Chart 
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Figure  44.  Speed  Advantage  Estimates  for  the  Loops 
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Figure  45,  Jennrich  Scramble  Microprogram  I 


il 
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Figure  46.  Jennrich  Scramble  M'croprogram  II 


Table  6 


FFT  TOTAL  TIME 


Array  Size 

Total  Time 

M 

N 

4 

16 

3  m  secs 

8 

256 

100  m  secs 

12 

4K 

2.3  sec 

16 

64K 

52  sec 

II. 5  ARRAY  OPTIMIZATION 


II.  5.1  Introduction 

An  aspect  of  array  design  based  upon  maximizing  gain  as  a  result  of  judicious 
choice  of  element  placement  in  accordance  with  a  model  of  the  noise  field  is  pre¬ 
sented.  The  technique  has  utility  in  a  relative  sense  since  its  practical  application 
is  dependent  upon  and  limited  to  situations  where  noise  data  is  available  and  a 
comparative  analysis  is  desired.  In  Section  2,  the  formulation  and  procedure  is 
outlined  and  some  results  are  presented  in  Section  3. 


II.  5. 2  Formulation 

Assuming  signal  coherence  and  equal  noise  power  at  each  element  in  an  N 
element  array,  the  conventional  processing  gain*’  can  be  expressed  as: 


G 


10  log 


N 


10 


1  +  (N-l)  p 


-dB. 


(83) 


where  the  average  correlation  coefficient  is: 


P 


N(N-l) 


N 

I 

i,J=l 


Pij 


(84) 


and  where  p^ 


is  the  pair-by-pair  element  noise  correlation  coefficient. 
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By  replacing  with  p(rj.),  data  derived  models  of  the  correlation  coeffi¬ 
cient  as  a  function  of  distance  Between  elements  (r^)  can  be  synthesized  for 
application  in  Equation  (83) .  By  restricting  the  class  of  functions  to  prohibit 
maximums,  the  average  value  of  the  differential  of  p  is: 


N 

d^  =(i/n)^s.  •  dr. , 

i=l 

where 


2 

N 

y 

~32<r  >~ 

r.. 

ij 

N(N-l) 

L 

j=i 

J 

L1  ij  U 

(85) 


(86) 


The  arrow  notation  is  used  to  denote  vector  quantities. 

A  2 

To  simplify  computation  p  has  been  modified  to  be  a  function  of  r..  as 

follows:  ^ 


S. 

i 


N  a,  2, 

4 _  V  8P(rij)  - 

N(N-l)  l  2  rij 

j=l  ij 

jA 


(87) 


This  eliminates  the  need  for  computing  |  7jj  |.  If  "ry  equals  zero  when  computing 
tq,  instrument  j  is  moved  to  reduce  the  magnitude  of  the  rj  which  is  measured 
relative  to  the  center  of  the  array.  This  restricts  two  instruments  from  having 
the  same  position. ^The  maximum  change  in  p  is  produced  when  drq  is  in  the 
same  direction  as  Sj;  therefore,  the  sensitivity  vector  Si  computed  for  each 
instrument  determines  the  direction  of  movement.  The  instrument  movement  is 
determined  by: 
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where  K  controls  the  incremental  step  size,  and,  thus,  the  magnitude  of  each 
instrument  move  is  proportional  to  the  ratio  of  its  sensitivity  vector  to  the  sys¬ 
tem  sensitivity. 


II. 5. 2.1  Constraints 

If  an  instrument  is  on  the  array  constraint  boundary,  an  instrument  movement 
outside  of  the  boundary  is  modified  to  be  the  projection  of  the  Aq  on  the  tangent 
to  the  constraint  circle  at  the  instrument  position.  When  an  instrument  is  moved 
outside  of  the  constraint  boundary,  the  instrument  is  returned  to  the  constraint 
boundary  along  the  instrument  vector  relative  to  the  center  of  the  array.  A  simi¬ 
lar  technique  is  used  for  other  types  of  constraints.  For  a  spoke  configuration 
where  all  instruments  must  lay  on  spokes  through  the  center  of  the  array,  in  addi¬ 
tion  to  being  within  a  fixed  diameter  circle,  the  instrument  position  vector  is 
rotated  to  correspond  to  the  nearest  spoke  after  each  move. 


II. 5. 2. 2  Program  Operation 

After  computing  the  directions  of  movement,  all  instruments  are  moved  at 
the  same  time.  Then,  a  new  set  of  sensitivity  vectors  is  computed.  This  process 
continues  until  either  the  sensitivity  vector  or  its  projection  on  the  constraint 
(when  constrained)  falls  below  a  variable  threshold.  A  three  instrument  configu¬ 
ration  is  shown  in  Figure  47  to  demonstrate  the  vector  relationship  associated 
with  the  program.  For  sufficiently  small  values  of  K,  this  technique  should  pro¬ 
duce  convergence  along  the  path  of  steepest  descent  toward  an  array  configuration 
with  a  minimum  value  of  p* . 


II. 5. 2. 3  Operator  Control 

To  permit  the  array  designer  to  dynamically  control  the  threshold  and  incre¬ 
ment  of  instrument  position  change,  the  array  configuration  is  displayed  on  an 
IBM  2250  CRT  Display.  Controls  for  modifying  the  program  parameters  are 
provided. 

The  array  configuration  that  satisfies  the  threshold  setting  is  a  function  of 
the  threshold  values.  A  better  instrument  configuration  can  be  achieved  by 
reducing  the  magnitudes  of  both  the  threshold  and  instrument  position  change. 

The  array  configuration  thus  achieved  satisfies  a  local  minimum  value  of  p*.  The 
program  does  not  have  a  test  to  determine  that  the  minimum  corresponds  to  a 
global  minimum  (the  ’’optimum"  configuration).  Confidence  in  selecting  an  array 
configuration  can  only  be  achieved  by  repeated  runs  (preferably  with  different 
initial  conditions)  to  test  for  better  design  configurations.  This  technique  is 
intended  only  to  provide  insight  into  selecting  an  array  geometry  by  providing  a 
model  to  serve  as  a  departure  point.  A  final  array  design  must  be  tempered  with 
geographic  constraints  and  data  acquisition  considerations. 
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II. 5. 3  Results 


II. 5. 3.1  Array  Performance 

Array  configurations  generated  using  the  optimization  program  are  shown 
in  Figure  48  for  circular  constraints  with  diameters  of  7  and  20  kilometers  using 
an  exponential  function  for  the  estimated  noise  correlation  model.  Predicted  array 
performance  using  exponential  functions  with  different  decay  rates  are  shown  in 
figures  49  and  50.  The  configurations  generated  by  the  array  optimization  pro¬ 
grams  are  compared  with  the  current  LASA  subarray  configurations  and  configu¬ 
rations  where  all  instruments  are  constrained  to  be  on  a  circle.  For  the  LASA 
seven-kilometer  subarray,  the  inner  rings  are  removed  to  demonstrate  the  change 
in  performance.  Using  both  noise  correlation  models,  the  LASA  E3  subarray  con¬ 
figuration  approached  the  maximum  gain  obtained  with  the  array  optimization  pro¬ 
gram.  The  LASA  subarray  configurations  are  shown  in  figures  51  and  52. 

The  16  kilometer  NORSAR  subarray  with  a  modified  hexagonal  geometry 
currently  being  installed  is  shown  in  Figure  53.  Perturbations  in  the  geometry 
resulted  from  practical  instrument  siting  considerations.  A  modified  circular 
geometry  and  the  configuration  generated  by  the  optimization  program  with  a 
16  kilometer  constraint  circle,  are  presented  in  figures  54  and  55,  respectively. 

The  modified  circular  geometry  also  resulted  from  practical  considerations.  In 
Figure  56,  the  optimal  gain  geometry  is  superimposed  on  the  NORSAR  subarray 
to  demonstrate  the  differences  in  the  configurations.  Table  7  summarizes  the  array 
performance  for  the  three  configurations  using  an  exponential  noise  correlation 
model.  The  NORSAR  subarray  and  the  circular  configuration  are  evaluated  over  a 
range  of  exponential  decay  rates.  In  addition,  the  minimum  communication  line 
length  required  in  each  of  the  configurations  is  tabulated. 


II. 5. 3. 2  Measured  Noise  Correlation  Data 

To  verify  the  validity  of  the  noise  correlation  models,  measurements  of  IASA 
subarray  performance  are  presented  in  figures  57,  58,  59,  and  60.  The  subarray 
performance  for  both  unfiltered  and  filtered  cases  is  shown.  These  performance 
characteristics  we  obtained  by  measuring  the  average  correlation  coefficient  for 
the  subarray  over  an  80  second  interval.  The  two  filters  used  were  third  order 
Butterworth  with  bandpasses  of  0.9  to  1.4  Hz  and  0.7  to  3.0  Hz,  respectively. 

Figure  61  shows  the  exponential  noise  correlation  models  used  for  evaluating  array 
performance.  These  models  can  be  compared  with  the  measured  correlation  coef¬ 
ficient  between  instruments  as  a  function  of  the  distance  between  them  for  three 
noise  records.  An  80  second  noise  record  was  used  for  each  case.  The  separate 
noise  records  are  identified  by  the  event  name  or  location  which  they  preceded. 

The  correlation  coefficients  were  fitted  with  a  sixth  order  polynomial  to  generate 
a  continuous  function.  This  data  is  presented  in  figures  62  through  70. 
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Figure  48.  Array  Configurations  for  20  Instruments 
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Figure  52.  LASA  Subarray 
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Figure  55.  Optimum  Gain  Geometry 
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Figure  56.  Superimposed  Geometry  (D  ~  16  Km) 
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Figure  57.  Observed  Subarray  Performance,  Montana  LASA  E3 
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Figure  58.  Observed  Subarray  Performance, 
Montana  LASA  I,  7  Km  Subarray 
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Figure  59.  Observed  Subarray  Performance, 
Montana  LASA  JJ,  7  Km  Subarray 
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Figure  60.  Observed  Subarray  Performance,  Montana  LASA  III,  7  Km  Subarray 
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Figure  61.  Noise  Correlation  Model 
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Figure  62.  Correlation  Coefficient  vs  Distance  Between  Instruments  I 


Noise  Prior  to  Kurile  Event  of  11/21/66 
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Figure  63  Correlation  Coefficient  vs  Distance  Between  Instruments  II 


Noise  Prior  to  Kurile  Event  of  11/21/66 
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Figure  64.  Correlation  Coefficient  vs  Distance  Between  instruments  III 


Noise  Prior  to  Kamchatka  Event  of  4/8/66 
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Figure  65.  Correlation  Coefficient  vs  Distance  Between  Instruments  IV 


Noise  Prior  to  Kamchatka  Event  of  4/8/66 
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Figure  66.  Correlation  Coefficient  vs  Distance  Between  Instruments  V 


Noise  Prio'  to  Kamchatka  Event  of  4/8/66 
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Figure  67.  Correlation  Coefficient  vs  Distance  Between  Instruments  VI 


Noise  Prior  fo  "Lonj'hot-"  Event  of  10^29/65 


00 

o 


>o  cn  o 

o  o  o  o 

luaioijjso^  uo|4D|9jjo^) 


CN 

o 

I 


123 


Figure  68.  Correlation  Coefficient  vs  Distance  Between  Instruments  VII 


Noise  Prior  to  "Longhot"  Event  of  10/29/65 
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Figure  69.  Correlation  Coefficient  vs  Distance  Between  Instruments  VIII 


Noise  Prior  to  "Longshot"  Event  of  10-29-65 
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Figure  70.  Correlation  Coefficient  vs  Distance  Between  Instruments  IX 


Appendix  III 


PHASE  DELAYS 


Operation  of  the  LASA  beam  former  requires  the  use  of  21  time  delays— one 
for  each  subarray.  Previous  studies  have  demonstrated  that  plane  wave  steering 
of  subarray  beams  must  be  augmented  by  correction  factors  to  take  account  of 
travel  time  anomalies. 

A  key  requirement  for  both  Detection  and  Event  Processing  is  the  library 
organization  of  phase  delays.  The  implications  of  using  a  library  form  of  process 
structure  are  examined  in  this  appendix.  Procedures,  based  on  correlating  a 
LASA  beam  with  each  of  the  subarray  beams  used  in  forming  it,  have  been 
developed  for  calculating  improved  steering  delays,  given  an  initial  set.  These 
procedures  described  in  Appendix  II— 1  will  be  employed  to  maintain,  correct,  and 
enlarge  the  library  of  phase  delay  parameters.  Regions  and  sectors  in  U-AZI-space 
are  defined,  and  associated  subarray  correction  factors  are  computed,  reflecting 
the  discrepancies  between  theoretical  and  actual  travel  times  observed  for  a 
sample  of  several  hundred  events.  A  test  performed  to  evaluate  the  effectiveness 
of  the  subarray  correction  factors  in  reducing  mis-steering  loss  is  also  described. 

The  methodology  to  cover  an  arbitrary  location,  U,  in  inverse  velocity  space 
by  a  judicious  choice  of  preformed,  presteered  subarray  beams  is  addressed  in 
Section  3.  Special  attention  is  given  to  subarray  beam  assignment  in  order  to 
minimize  the  expected  loss. 


III.l  PROCESS  STRUCTURE 

The  Detection  Processor  forms  LASA  beams  by  summing  in  phase  appropriate 
beams  chosen  from  a  set  of  preformed,  pre- steered  subarray  beams.  Since  the 
delays  required  for  steering  each  subarray  are  a  function  of  the  location  of  the 
LASA  beam  being  formed,  they  will  be  obtained  from  a  library  of  phase  delays, 
organized  by  LASA  beam  coordinates.  In  Event  Processing  the  library  delays 
will  be  employed  or  the  optimal  delays  will  be  obtained  from  correlations  that  can 
be  generated  for  events  of  special  interest. 
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III.  1.1  Phase  Delay  Library  Parameters 

The  entries  for  the  phase  delay  library  parameters  will  be  a  combination  of 
several  source  factors  among  which  are 

a.  Regular  plane  wave  steering  delays 

b.  Sets  of  subarray  steering  corrections  each  of  which  is  valid  over  a 
specified  region  of  U-AZI  space 

c.  Inter-regional  corrections  calculated  from  the  sets  of  subarray 
corrections 

d.  Optimal  delays  obtained  as  a  result  of  an  iterative  process  based  on 
correlation  methods. 

A  library  structure  demands  that  requirements  be  specified  for  establishing, 
maintaining,  and  updating  library  entries. 

It  is  planned  to  employ  plane  wave  delays  for  beam  steering,  particularly  for 
steering  to  aseismic  locations. 

For  seismic  areas,  it  is  planned  to  modify  the  plane  wave  delays  by  means 
of  empirically  derived  subarray  correction  factors;  each  set  is  valid  over  a  speci¬ 
fied  region  in  U-AZI  space.  The  development  of  35  sets  of  regional  correction 
factors  will  be  described  in  this  appendix. 

Formulas  have  been  developed  for  interpolating  between  regions  in  order  to 
obtain  correction  factors  for  events  located  outside  regions.  Plane  wave  delays 
will  be  replaced  by  optimal  delays  for  certain  regions  and  for  weak  events.  These 
will  be  obtained  by  use  of  the  correlation  program  described  in  Appendix  II- 1. 

It  will  be  routine  practice  to  check  the  subarray  regional  corrections  and  the 
interpolation  formulas  on  a  regular  basis.  This  will  detect  any  systematic  devia¬ 
tions  which  may  reflect  trends  in  the  correction  factors. 

Analysis  of  strong  events  will  be  the  chief  instrument  used  to  correct  and 
enlarge  library  phase  delays.  Detection  and  identification  of  a  strong  event  will 
be  followed  by  the  input  of  its  delays  to  the  correlation  program,  which  will  calcu¬ 
late  optimal  delays  for  the  event.  The  optimal  delays  may  be  used  to  check 
existing  region  correction  factors  and  associated  interpolation  formulas,  or  may 
serve  as  the  initial  set  of  correction  factors  for  a  new  region. 

Criteria  must  be  established  for  determining  when  the  correlation  program 
is  to  be  used  for  updating  library  entries,  and  for  determining  when  an  event  may 
be  considered  so  strong  that  its  delays  may  be  used  to  correct  library  parameters. 
The  interplay  between  library  parameters  and  the  correlation  program  is  shown 
in  Figure  71.  For  routine  events,  appropriate  phase  delays  will  be  obtained  from 
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Figure  71  .  Interaction  between  Phase  Delays  and  the 
Correlation  Program 
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the  library.  For  special  events,  in  particular,  weak  ones,  the  correlation  program 
will  be  employed  to  calculate  optimal  delays,  using  as  input  either  modified  plane 
wave  delays  or  delays  for  the  beam  on  which  the  event  was  detected. 


III. 2  ANOMALY  CHARACTERIZATION 


III.  2.1  Introduction 

This  section  describes  the  development  of  steering  correction  factors 
resulting  from  an  examination  of  several  hundred  events,  each  of  which  yielded  a 
set  of  21  time  delay  corrections.  It  is  shown  that  these  sets  of  correction  factors 
can  be  grouped  to  produce  average  sets  of  corrections  which  are  valid  over  judi¬ 
ciously  selected  geographic  regions.  Because  of  the  irregular  space  and  time 
distribution  of  earthquakes,  these  regions  do  not  cover  the  entire  earth.  Conse¬ 
quently,  it  is  necessary  to  interpolate  to  obtain  correction  factors  for  beams  that 
are  being  steered  to  areas  located  between  the  nominally  selected  regions.  The 
development  of  such  interpolation  formulas  is  described  in  this  appendix. 


III. 2. 2  Method  of  Correction 

Theoretical  steering  delays  must  be  amended  by  adding  some  correction  fac¬ 
tor.  This  is  currently  being  accomplished  by  defining  for  each  event  region  and 
receiver  site  a  travel  time  anomaly.  This  anomaly  is  the  difference  between  the 
actual  travel  time  and  a  theoretical  travel  time,  based  on  a  worldwide  travel  time 
average  for  events  at  a  given  range,  such  as  those  found  in  the  Jeffreys  Bullen 
Seismological  Tables.  The  required  event  range  can  be  obtained  for  each  event 
from  an  event  location  based  on  the  world  seismic  network.  Time  of  origin  is 
also  determined  by  this  network.  The  average  anomaly  for  each  receiver  for  a 
given  geographic  region  is  then  calculated  from  a  number  of  events  occurring  in 
that  region,  and  is  used  to  provide  a  correction  to  the  theoretical  travel  time 
tables  to  give  the  required  beam  steering  delays.  The  LASA  system  steering 
error  for  a  given  event  would  depend  on  how  well  the  event  arrival  time  differences 
fit  the  region  averages.  An  interesting  and  informative  study  of  time  anomalies 
using  this  method  on  the  LASA  array  has  been  completed  by  E.  F.  Chiburis.^^ 

Another  possible  method  of  organizing  the  time  delays  is  to  establish  a  best 
fitting  plane  or  second  order  wavefront  for  each  event,  and  to  use  as  anomalies 
the  variabilities  of  the  deviations  from  this  best  fitting  wavefront.  As  in  the 
currently  used  method,  the  delays  are  based  on  the  anomaly  averages  for  a  given 
region. 

A  program  to  implement  the  wavefront  method  has  been  written.  This  pro¬ 
gram,  called  the  Least  Squares  Plane  and  Quadratic  Wavefront  Program,  accepts 
as  inputs  the  positions,  arrival  times,  and  computation  weights  of  up  to  25 
seismometers  (or  subarray  beams),  and  calculates  three  different  least  squares 
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wavefronts  with  the  actual  arrival  time  deviations  from  the  best  fits.  The 
three  fits  are 

a.  A  best-fit  plane  wavefront 

b.  A  best-fit  plane  wavefront  followed  by  an  adjustment  of  arrival 
times  to  consider  the  expected  effect  of  range  upon  wavefront 
curvature,  with  the  above  process  repeated  for  any  desired 
number  of  iterations 

c.  A  best-fit  quadratic  wavefront  plus  the  three  derived  second 
derivatives  or  curvatures,  expressed  in  array  coordinates, 
range  coordinates,  and  the  principal  coordinates  of  the  quadratic 
surface. 

The  input  data  for  the  wavefront  program  consisted  of  optically  read  arrival 
time  data,  made  available  by  Seismic  Data  Laboratories  (SDL)  of  Teledyne  Indus¬ 
tries,  Inc.,  and  the  VELA  Seismological  Center  (Air  Force  Technical  Application 
Center),  Alexandria,  Virginia.  A  total  of  approximately  four  hundred  event 
arrivals  were  made  available  in  their  original  form  by  SDL.l 7 


ni.2.3  Procedure  for  Obtaining  Steering  Delay  Corrections 

Chiburis  describes  the  grouping  of  several  hundred  events  into  36  geographic 
regions.  Table  8,  which  uses  the  regions  so  defined  together  with  later  data  from 
SDL,  lists  the  number  of  events  in  each  region.  As  described  in  Reference  16,  a 
region  is  formed  by  grouping  teleseisms  whose  subarray  travel  time  anomalies 
are  consistent  within  the  region.  The  377  events  and  the  36  regions  form  the 
starting  point  of  the  development  reported  herein. 

The  377  events  have  been  processed  through  the  Least  Squares  Wavefront 
Program  and  the  resulting  computed  values  of  azimuth  (AZI)  and  inverse  phase 
speed  (U)  have  been  plotted  for  each  event  as  shown  in  Figure  72.  The  other  out¬ 
put  of  the  wavefront  program— the  deviation  in  seconds  of  the  actual  arrival  from 
the  best  fitting  plane  wavefront— serves  as  input  to  the  Seismic  Steering  Delay 
Anomalies  Program.  Given  a  group  of  events,  each  with  its  set  of  21  correction 
factors  or  deviations  (one  for  each  subarray),  this  program  calculates  the  average 
deviation  at  each  subarray.  The  difference  between  each  eventfs  deviation  and  the 
average  deviation  of  the  group  at  each  subarray  is  computed,  yielding  a  set  of  21 
differences,  i.e.,  deviations  of  deviations.  The  average  and  standard  deviation  of 
this  set  of  differences  and  the  mis-steering  loss  in  dB  for  the  event's  wavefront 
are  calculated. 

The  division  of  the  377  events  into  36  regions^  furnished  the  initial  groups 
used  in  the  Steering  Delay  Anomalies  Program.  The  program  results  indicated 
that  some  of  the  events  yielded  excessive  dB  losses  and  it  was  decided  to  regroup 
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Table  8 


377  EVENTS  GROUPED  INTO  36  REGIONS 


Region 

No. 

Geographic  No.  of  Events 

Area  In  Region 

Region 

No. 

Geographic  No.  of  Events 

Area  in  Region 

1 

Central- So. 
Alaska 

6 

19 

N.  Atlantic 

2 

2 

Kodiak  Island 

11 

20 

Azores 

2 

3 

Alaska  Pen. 

3 

21 

No.  West  Indies 

5 

4 

Unimak 

4 

22 

Venezuela- So. 

West  Indies 

3 

5 

Fox- Andreanof 

19 

23 

Eastern  Mexico 

9 

6 

Andreanof 

4 

24 

N.  Centr.  Ameria 

8 

7 

Andr-Rat 

22 

25 

S.  Centr.  America 

8 

Near  Island 

14 

26 

N.  Columbia 

12 

9 

Kamchatka 

12 

27 

Peru-Brazil 

8 

10 

Kurile 

29 

28 

Peru-Bolivia 

13 

11 

Sakhal  in- 
Hokkaido 

8 

29 

N.  Chile- 
Bolivia 

16 

12 

Honshu 

Marianas 

20 

30 

N.  Chile- 
Argentina 

12 

13 

Marianas- 

Carolines 

13 

31 

C.  Chile- 
Argentina 

13 

14 

Taiwan-Ryuku 

10 

32 

West  Mexico 

4 

15 

Kazakh 

9 

33 

Easter  Island 

3 

16 

Hindu  Kush 

7 

34 

Easter  Island 

7 

17 

East  Caucasus 

3 

35 

Samoa- Tonga 

22 

18 

Greece -Turkey 

14 

36 

Fiji  Islands 

23 
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Figore  72.  3 77  Wavefronts,  Calculated  U  and  AZi  as  Observed 
at  LASA 
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the  events  into  new  regions.  There  are  35  regions  in  the  new  grouping.  They 
are  plotted  in  Figure  73  and  their  coordinates  are  listed  in  Table  9.  The  revised 
regions  are  quite  similar  to  the  original  36.  For  example,  revised  regions  109 
through  123  correspond  exactly  to  original  regions  12  through  26  respectively. 

After  regrouping  the  377  events  into  35  new  regions,  averages  were  again 
formed  and  the  anomalies  program  was  used  to  check  the  assignment  of  the 
events  to  the  regions  and  to  assign  weights  to  the  events.  For  those  events 
which  continued  to  show  a  loss  of  more  than  1  dB  at  1  Hz,  a  weight  of  zero  was 
assigned  while  the  others  were  weighted  one.  New  averages  were  formed  for 
each  region,  using  only  those  events  weighted  one,  i.e. ,  using  those  events  whose 
loss  was  less  thanldB  at  1  Hz.  The  anomalies  program  was  used  to  provide  a 
final  check.  This  time  there  were  no  further  changes— the  events  weighted  one 
continued  to  show  a  loss  of  less  than  1  dB  at  1  Hz  while  those  weighted  zero 
showed  even  greater  losses  (as  was  to  be  expected  since  they  had  been  excluded 
from  the  average  forming  prooess  the  second  time).  Thus,  the  phase  delay  aver¬ 
ages  for  the  revised  regions  are  arrived  at  by  a  two  step  iterative  process,  based 
on  stable,  similar  events  each  of  which  experiences  a  loss  of  less  than  1  dB  at  1  Hz 
due  to  steering  to  these  regional  averages.  Table  10  lists  the  number  of  events 
with  weight  one  for  each  region.  As  shown  in  the  table,  there  is  a  total  of  335 
weight  one  events  out  of  a  total  of  377  events  and  only  42  events  were  weighted 
zero. 

The  Anomalies  Program  yields  the  average  loss  (in dB)  per  region,  i.e.,  the 
standard  deviations  of  the  events  in  a  region  are  averaged,  and  converted  to  a 
loss. 


For  an  event  occurring  in  a  region,  this  average  loss  figure  represents  the 
loss  in  dB  which  results  from  using  the  regional  averages  as  steering  correction 
factors  rather  than  a  set  of  corrections  that  would  be  optimum  for  the  particular 
event.  Table  11  lists  the  average  loss  per  region  for  the  35  revised  regions.  As 
can  be  seen  from  the  table,  che  average  loss  is  less  than  1  dB  at  1  Hz  for  every 
region.  These  results  indicate  that  the  regional  averages  can  be  instrumental  in 
providing  correction  factors  for  beams  deployed  within  regions. 


in. 2. 4  Definition  of  Sectors 

To  obtain  correction  factors  for  areas  having  no  history  of  teleseismic 
activity,  some  form  of  interpolation  between  regions  must  be  developed.  Inter¬ 
polation  between  regions  is  based  on  the  assumption  that  correction  factors  are 
functions  only  of  azimuth  or  inverse  phase  speed  or  both.  This  implication  or 
assumption  must  be  considered  somewhat  questionable,  for  examination  of  the 
regional  corrections  reveals  no  discernible  pattern  in  the  average  deviations. 
They  appear  not  to  be  simple  functions  of  azimuth  or  of  range,  but  they  fluctuate 
in  an  erratic  manner.  Figures  74,  75,  and  76  illustrate  the  deviations  per  region 
for  three  specific  subarrays. 
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118 


egic 

No. 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 


Table  10 


NUMBER  OF  CONTRIBUTING  EVENTS 


No.  of 
Events 

6 

10 

7 
37 
13 
16 
13 
16 
19 
11 

8 
9 


Region 

No. 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 


No.  of 
Events 

5 

2 

12 

2 

2 

5 

3 
9 
8 
7 

10 

4 


Region 

No. 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 


No.  of 
Events 

8 

6 

8 

13 

10 

7 

4 

5 
2 
4 

34 
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Table  11 


AVERAGE  LOSS  PER  REGION 


Region 

No.  of  Events 

Weight  One /Total 

Average  Loss 
(in  dB) 

101 

6/6 

0.44 

102 

10/13 

0.55 

103 

7/7 

0.48 

104 

37/41 

0.45 

105 

13/16 

0.34 

106 

16/17 

0.39 

107 

13/14 

0.33 

108 

16/18 

0.25 

109 

19/20 

0.36 

110 

11/13 

0.34 

111 

8/10 

0.31 

112 

9/9 

0.39 

113 

5/7 

0.49 

1 14 

2/3 

0.43 

115 

12/14 

0.45 

116 

2/2 

0.15 

117 

2/2 

0.15 

118 

5/5 

0.36 

119 

3/3 

0.77 

120 

9/9 

0.53 
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Table  11 


Region 


121 

122 

123 

124 

125 
12G 

127 

128 

129 

130 

131 

132 

133 

134 

135 

Total 


(concluded) 

No.  of  Events 
Weight  One/Total 

8/8 

7/7 

10/12 

4/4 

8/9 

6/7 

8/8 

13/14 

10/13 

7/7 

4/4 

5/7 

2/3 

4/5 

34/40 

335/337 


Average  Loss 
(in  (IB) 


0.45 

0.27 


0.38 


0.36 


0.41 


0.36 

0.32 


0.39 


0.44 

0.41 

0.21 


0.54 


0.16 

0.41 


0.36 


Figure  74  Regional  Deviations  (seconds)  for  Subarray  Cl 
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Figure  75.  Regional  Deviations  (seconds)  for  Subarray  C2 
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An  interpolation  procedure  has  been  developed  in  an  attempt  to  estimate 
correction  factors  for  areas  near  regions.  The  interpolation  method  is  based 
on  the  construction  of  areas  called  sectors  which  are  essentially  defined  in  terms 
of  regions.  A  sector  is  an  area  containing  two  or  more  regions.  The  guideline 
used  in  establishing  sectors  has  been  to  define  their  boundaries  so  that  a  sector 
is  essentially  constant  either  for  azimuth  or  inverse  phase  speed.  Figure  77 
depicts  the  nine  sectors  that  have  been  defined. 

Table  12  gives  the  U-AZI  coordinates  tor  each  sector  and  lists  the  regions 
contained  in  each  sector.  As  shown  in  Figure  78  sectors  do  not  extend  significantly 
beyond  the  regions  they  contain. 


III. 2. 5  Methods  of  Interpolation  and  Extrapolation 

The  interpolatory  method  used  to  obtain  a  deviation  for  each  subarray  and 
sector  is  simple  linear  interpolation  in  which  either  azimuth  or  speed  is  the 
independent  variable.  For  example,  if  an  event  occurs  in  sector  A,  the  deviation 
to  be  used  for  steering  to  that  event  is  obtained  by  interpolating  in  azimuth  between 
regions  114  and  115.  The  linear  interpolation  formula  is  stated  in  Table  13. 

Table  14  gives  the  (X,  Y)  values  relative  to  subarray  B1  for  each  region  in  each 
sector.  Here,  Y  is  the  value  of  the  deviation  while  X  is  either  an  azimuth  value  or 
a  U  value,  depending  on  whether  the  independent  variable  is  azimuth  or  inverse 
phase  speed.  The  tables  for  the  other  subarrays  are  similar  and  therefore  not 
presented. 

Two  or  more  regions  are  sometimes  combined  and  treated  as  a  single  region. 
In  particular,  regions  122  and  124,  which  are  contained  in  Sector  C,  are  assigned 
the  common  azimuth  value  of  147°  and  the  common  deviation  of  0.017  seconds;  this 
pair  of  values  is  used  for  all  21  subarrays. 

The  interpolation  procedure  will  be  demonstrated  by  computing  the  correction 
factor  for  subarray  B1  steered  to  event  V  whose  location  is  estimated  at  U  =  0.080 
and  AZI  =  130°.  Figure  78  shows  that  V  is  in  sector  C  but  not  within  any  region 
in  C.  Table  14  shows  that  deviation  is  a  function  of  azimuth  in  sector  C;  with 
respect  to  that  variable,  V  is  located  between  regions  119  and  123.  Employing 
the  notation  of  Table  13,  let  (X^,  Y^)  represent  the  coordinates  of  region  119  and 
(X2,  Y2)  represent  the  coordinates  of  region  123. 

Thus,  (Xr  Yx)  =  (123°,  -  0.062), 

(X2,  Y2)  =  (137°,  -  0.188) 

and  (X,  Y)  =  (130°,  Y),  since  AZI  is  the  independent  variable  in  this  case.  The 
deviation  Y  is  computed  as  follows: 
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Table  12 


SECTORS 


Sector  Range  Window  Geographic  Azimuth 

(sec/Km)  Window  (degrees) 


A 

0.035  < 

U< 

0.052 

320  <  0  <  45 
(Through  360) 

B 

0.060  < 

U< 

0.070 

45  <  0<  75 

C 

0.067  < 

U  < 

0.090 

110  <  0  <  170 

D 

0.036  < 

U  C 

0.067 

140  <  0  <  170 

E 

0.036  < 

U  < 

0.090 

170  <  0  <  190 

F 

0.036  < 

U  < 

0.048 

225  <  0  <  250 

G1 

0.039  < 

U  < 

0.090 

290  <  0  <  304 

G2 

0.039  < 

U  < 

0.090 

304  <  0  <  312 

G3 

0.039  < 

U  < 

0.090 

312  <  0  <  320 

Regions  Contained 
in  Sector 

112,  113,  114,  115 
116,  117 

118,  119,  120,  121 
122,  123,  124 

125,  126,  127,  128 
129,  130 

131,  132,  133 

134,  135 

103,  104,  110 

102,  105,  109 

101,  106,  107,  108 
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Figure  78,  Relationship  Between  Regions  and  Sectors 
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Table  13 


LINEAR  INTERPOLATION  FORMULA 


Given:  Points  (Xr  and  (X„,  Y.,) 


The  straight  line  Y  =  aX  +  b  which  passes  through  (X^ ,  Y 
may  be  expressed  as: 


Y  =  Y  + 


X  -  x1 

X2~X1 


(Y2  -  Yl> 


j)  and  (X2,  Y2) 
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Table  14 


SUBARRAY  FORMULAS  FOR  SECTORS 
SUBARRAY  B1 


Sector 

Region 

U 

(Sec/Km) 

Azimuth 

(degrees) 

Deviation 

(seconds) 

A 

112 

0.0435 

348 

0.144 

113 

0.0435 

355 

0.147 

114 

0.0435 

22 

0.210 

115 

0.0435 

39 

0.147 

B 

116 

0.065 

52 

0.090 

117 

0.065 

73 

0.073 

C 

118 

0.078 

119 

-0.011 

119 

0.078 

123 

-0.062 

123 

0.078 

137 

-0.188 

122-124 

0.078 

147 

-0.209 

121 

0.078 

154 

-0.253 

120 

0.078 

161 

-0.273 

D 

125 

0.064 

145 

0.238 

126 

0.060 

145 

0.184 

127 

0.057 

145 

0.131 

128 

0.054 

145 

0.078 

129 

0.050 

145 

0.024 

130 

0.046 

145 

0.066 

E 

131 

0.085 

185 

-0.159 

132 

0.061 

185 

-0.075 

133 

0.054 

185 

-0.012 

F 

134 

0.042 

224 

0.151 

135 

0.042 

224 

0.151 

G1 

103 

0.077 

297 

0.054 

104 

0.070 

297 

0.098 

110 

0.043 

297 

0.044 

G2 

102 

0.082 

308 

0.041 

105 

0.067 

308 

0.031 

109 

0.048 

308 

0.016 

G3 

101 

0.081 

316 

0.100 

106 

0.064 

316 

0.001 

107 

0.058 

316 

-0.053 

108 

0.054 

316 

0.011 

147 


(X  -  xx) 


Y  Y1  +  (X2  -  Xx)  '  (Y2  "  Yl* 


Y  =  -0.062  +  (-0-188  -  (-0.062))  =  -0.125  seconds. 


III.2.6  Uncovered  Areas 

It  is  also  necessary  to  determine  deviations  for  events  which  occur  in  areas 
not  contained  in  any  sector.  Since  such  areas  are  not  near  areas  in  which  numerous 
events  are  observed,  little  knowledge  of  the  expected  phase  delays  exists.  One 
possible  solution  is  to  steer  by  means  of  plane  waves  with  no  subarray  deviation 
correction.  A  second  possible  method  is  to  interpolate  or  extrapolate  using  sets 
of  regional  averages  from  several  regions.  A  third  possible  solution  is  to  define 
more  complex  functions  which  will  cover  the  entire  teleseismic  zone.  This  method 
appears  most  promising,  but  is  contingent  upon  the  results  achieved  by  the  pres¬ 
ently  defined  nine  sectors. 


III. 2. 7  Testing  The  Regions  And  Sectors 

This  section  discusses  a  test  which  was  performed  to  evaluate  steering  delay 
correction  factors  based  on  the  construction  of  regions  and  sectors.  The  results 
of  the  test  indicate  that  regional  correction  factors  can  be  employed  successfully 
to  obtain  improved  phase  delays. 

The  test  commences  with  a  simulated  steering  of  the  array  to  an  event  located 
within  a  region,  by  means  of  plane  wave  delays  modified  by  correction  factors 
associated  with  the  region.  Next,  a  best  fit  plane  wave  is  obtained  for  the  event 
by  using  the  time  of  arrival  at  the  center  seismometer  of  each  subarray  as  input 
data  for  the  Least  Squares  Wavefront  Program.  The  difference  between  computed 
and  actual  arrival  time  at  the  center  seismometer  for  each  subarray  is  next  cal¬ 
culated.  Finally  the  Seismic  Steering  Delay  Anomalies  Program  is  used  to  evalu¬ 
ate  the  difference  (if  any)  between  the  steering  correction  factor  and  the  best  fit 
plane  wave  deviation  at  each  subarray.  The  Anomalies  Program  calculates  from 
the  21  differences  the  loss  (in  dB)  that  can  be  attributed  to  using  subarray  steering 
correction  factors.  The  magnitude  of  this  loss  measures  the  effectiveness  of  the 
set  of  steering  corrections  associated  with  the  region  containing  the  event.  Since 
a  region  is  characterized  by  its  set  of  steering  corrections,  this  loss  also  gauges 
the  reliability  of  the  region  itself.  If  the  loss  is  unacceptably  large  (where  unac¬ 
ceptably  large  would  have  been  previously  specified),  the  conclusion  is  that  one  or 
more  of  the  following  factors  is  responsible. 

a.  The  region  is  not  well-defined,  i.e. ,  the  correction  factors  are 
ineffective  for  steering  the  array. 

b.  The  event  was  inaccurately  timed  at  several  subarrays . 
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c.  A  grossly  inaccurate  reading  was  made  at  one  or  more  subarrays. 

d.  Some  other  parameter,  such  as  depth,  must  be  taken  into  account. 

Factors  b,  c,  and  d  can  be  eliminated  by  rechecking  seismometer  records  and 
repeating  the  computations  as  required.  If  a  recalculation  continues  to  yield 
an  unacceptably  large  loss,  it  must  be  concluded  that  the  steering  correction 
factors  are  of  no  benefit  in  reducing  mis-steering  loss.  Such  a  test  procedure 
will  be  more  meaningful  and  reliable  if  there  are  several  events  per  region  which 
can  be  tested. 


III. 2. 8  Test  Results 

The  test  conducted  in  this  quarter  was  based  on  a  set  of  56  events,  each  with 
center  seismometer  arrival  times,  made  available  by  Seismic  Data  Laboratory, 
Teledyne, Inc. 

The  events  were  processed  through  the  Least  Squares  Wave  Front  Program 
which  calculated  inverse  phase  speed  and  direction  for  each  event.  As  shown  in 
Figure  79,  it  was  found  that  33  events  fell  inside  regions  and  23  fell  outside. 

The  results  of  the  test  for  the  33  events  which  fell  inside  regions  are  shown 
in  Table  15.  Columns  1  and  2  are  self-explanatory.  Column  3  indicates  the  num¬ 
ber  of  center  seismometers  that  were  used  in  the  curve  fitting  process.  The 
fewer  the  number  used,  che  less  reliable  is  the  fit.  Inspection  of  the  loss  figures 
given  in  Column  4  shows  that  26  of  the  33  events  suffered  a  loss  of  1  dB  or  less 
at  1  Hz.  Thus  it  may  be  concluded  that  regional  correction  factors  would  have 
helped  steer  to  these  events  successfully.  As  for  the  remaining  7,  it  is  seen 
from  Column  3  that  for  4  of  them,  the  number  of  sensors  used  was  10  or  less  out 
of  21.  This  small  number  of  sensors  implies  that  the  wave  fitting  is  not  to  be 
regarded  as  valid  or  reliable  because  the  deviations  from  the  best-fit  wave  front 
are  so  large  that  they  cannot  be  compared  with  regional  averages.  Consequently 
these  events  have  been  disregarded. 

Interpolation  procedures  were  not  available  for  application  to  the  23  events 
falling  outside  the  regions.  However,  some  of  these  events  lay  rather  close  to 
regions  and  provided  an  opportunity  to  examine  how  far  regional  deviations  could 
be  extended.  Accordingly,  it  was  decided  to  simulate  steering  to  eight  of  these 
events  using  for  the  simulation  the  correction  factors  for  the  regions  closest 
to  the  events.  These  eight  events  are  shown  in  Figure  80;  the  closest  regions  are 
shaded  in  the  figure.  The  results  of  this  attempt  to  stretch  regional  corrections 
are  shown  in  Table  16.  Columns  1  through  4  are  similar  to  Columns  1  through  4 
of  Table  16.  Column  5  contains  relevant  comments.  As  can  be  seen  from  Table  16, 
five  of  the  eight  could  have  been  steered  to  successfully  using  regional  averages. 

The  three  events  for  which  regional  average  steering  would  have  been  unsuccessful 
have  been  labeled  by  a  check  (\J)  in  Figure  80. 
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Figure  79.  23  Events  Located  Outside  the  35  Regions 
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Table  15 


TESTING  OF  THE  REGIONAL  AVERAGES  USING 
33  EVENTS  LOCATED  WITHIN  REGIONS 


Region 

Number 

Event 

Identification 

Number 

Number 

of 

Sensors 

Loss 
(in  dB) 

101 

4840 

6 

1.65 

4850 

13 

1.41 

4880 

20 

0.35 

5000 

20 

0.55 

5020 

18 

0.37 

5050 

18 

0.54 

103 

4830 

17 

0.33 

4980 

17 

0.60 

5100 

20 

0.57 

104 

4635 

20 

0.87 

4655 

20 

0.96 

4680 

20 

0.20 

5145 

19 

0.30 

5155 

18 

0.35 

106 

4625 

20 

0.51 

108 

4310 

19 

0.23 

110 

4639 

9 

1.0! 

114 

4930 

20 

0.70 

115 

4780 

16 

1.01 

4810 

14 

0.43 

4860 

17 

0. 57 

4940 

21 

0.50 

116 

4627 

10 

1.41 

121 

4700 

11 

0.74 

4790 

19 

0.81 

4870 

10 

1.31 

124 

4920 

21 

0.26 

4950 

21 

0.32 
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Table  15 
(Concluded) 


Region 

Number 

Event 

Identification 

Number 

Number 

of 

Sensors 

Loss 
(in  dB) 

125 

5060 

21 

0.37 

128 

4630 

21 

0.23 

131 

4800 

19 

0.34 

5010 

12 

1.45 

135 

4320 

20 

0.36 
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Figure  80„  Steering  to  8  Events  Near  Regions 
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Table  1G 


TEST  OF  STEERING  TO  EVENTS  LOCATED  NEAR 


REGIONS  USING  REGIONAL  AVERAGES 

Region 

Event 

Number 

Loss 

Comment 

Number 

Identification 

of 

(in  dR) 

Number 

Sensors 

101 

4960 

21 

1.38 

102 

4623* 

11 

1.61 

This  set 
unacceptable 

105 

4623* 

11 

0.71 

This  steering 
acceptable 

109 

4629 

13 

0.23 

5115 

20 

0.23 

116 

4631 

17 

0.52 

118 

5080* 

18 

i.  ig'I 

Neither  set  is 
acceptable 

119 

5080* 

18 

1.73  f 

129 

5070* 

18 

1.01 

Unacceptable 

130 

5070* 

18 

0.79 

Acceptable 

132 

4637 

15 

3.24 

*Event  is  steered  to  using  two  or  more  sets  of  regional  averages. 
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III. 2. 9  Conclusions 


The  test  results  verify  that  the  present  regions  are  well-defined.  Correction 
factors  associated  with  the  regions  have  been  used  to  steer  successfully  to  nearlv 
all  events  located  within  the  regions.  The  fact  that  five  events  located  near 
regions  were  also  steered  successfully  by  using  the  regional  corrections  lends 
credence  to  the  planned  method  of  interpolating  between  regions. 
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III.  3  SUB  ARRAY  ASSIGNMENT 


III.  3.1  Introduction 

The  task  addressed  here  is  to  cover  an  arbitrary  location,  Uj  ,  in  inverse 
velocity  space  by  a  judicious  choice  of  pre-formed,  pre-steered  subarray  beams. 
Ordinarily  the  problem  is  resolved  by  choosing  that  set  of  beams,  one  from  each 
subarray,  which  is  closest  to  the  point  Uj.  However,  consider  the  events  plotted 
in  Figures  81,  82,  and  83.  The  square  (  □)  in  each  figure  denotes  the  location  of 
the  event  in  question.  Observed  seismometer  arrival  times  were  used  to  best-fit 
a  plane  wave  to  each  of  the  21  subarrays.  ^Then  for  each  subarray,  the  best-fitting 
plane  wave  was  used  to  calculate  a  value  U^j  which  may  be  interpreted  as  the 
event  location  as  estimated  by  using  only  subarray i.  For  each  event,  the  21 
values  Ugj  have  also  been  plotted  in  Figures  81,  82,  and  83_where  they  are  shown 
as  dots  ( . ).  It  is  seen  from  the  figures  that  in  no  case  is  equal  to  U-  for  any 
pair  (j*£). 

The  situation  depicted  in  these  figures  has  several  implications  for  beam¬ 
forming  and  subarray  steering.  Since  the  dots  denote  the  positions  to  which  the 
subarrays  should  be  steered  to  form  the  best  LASA  beam  covering  the  correspond¬ 
ing  square  (□  ),  it  is  possible  that  for  one  or  more  subarrays  the  preformed  beam 
which  is  closest  to  the  nominal  location  Uj  is  not  the  same  as  the  beam  which  is 
closest  to  U^j.  If  this  is  true,  the  assignment  of  subarray  beams  is  no  longer  an 
elementary  task. 

The  fact  that  Uj  and  U^j  may  not  be  the  same  point  indicates  that  in  some 
cases  phase  delays  will  be  obtained  by  a  chaining  process.  The  location  being 
steered  to  will  not  supply  the  delays,  but  it  will  point  to  several  other  locations, 
each  of  which  will  furnish  its  associated  delay.  Organization  of  the  phase  delay 
library  must  take  cognizance  of  this  possibility. 


III. 3.2  Notation  and  Definitions 

Notation  for  indexes, as  used  in  this  section,  is  defined  as  follows: 


NOTATION 

DEFINITION 

RANGE 

i 

instrument  index 

1 ,  •  .  •  ,  I 

j 

LASA  arr&y  beam  index 

1 ,  •  •  •  ,  J 

k 

subarray  beam  index 

1,  .  .  .  ,  K 

i 

subarray  index 

1,  ....  E 
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Figure  81  „  Kamchatka  Event 
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Figure  82„  Kazakh  Event 
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Figure  83  o  Novaya  Event 
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The  following  bounds  may  be  assumed  for  the  above  indexes:  I  £  25,  400  £  J  £  1024, 
K  £  10,  and  L  <  21.  The  inequalities  I  <  25  and  L  <  21  imply  that  selective  assign¬ 
ment  of  instruments  and  subarrays  may  be  utilized  for  beamforming  processes. 


Predefined  sets  of  points  in  inverse  velocity  space,  assumed  as  inputs  for 
the  problem  of  subarray  assignment,  are  as  follows: 


a.  The  arbitrary  set  {Uj  =  (Uxj ,  Uyj)  }  ,  assumed  to  be  covered 

by  a  corresponding  set  of  LASA  beams. 


b.  The  set  {  U  =  (UXj^,  Uy^j)}  ,  where  U^j  is  interpreted  as  an 
"aiming  point"  covered  by  that  beam  of  subarray  i  which  best 
contributes  to  the  formation  of  LASA  beam  j,  under  optimal 
processing.  This  set  of  points  is  referenced  to  the  set{  Uj}  , 
each  point  being  an  estimated  location  of  TTj ,  as  viewed  from 
subarray  f.  It  is  noted  that  in  the  absence  of  anomali 


sub array 
for  each  pair  i,  j. 


alies,  U.:  =  TT. . 
J  ij 


c.  The  set  {Ufk  =  U  ik)},  where  Ufk  is  interpreted  as  an 

"aiming  point"  covered  Dy  the  K^h  beam  of  subarray i.  This  set 
of  points  is  assumed  to  have  been  selected  on  the  basis  of  some 
previously  defined  optimization  criteria  and  is  subject  to  reassign¬ 
ment  on  the  basis  of  these  criteria. 


III. 3. 3  Evaluation  of  Loss 


Evaluating  the  maximum  loss  observed  for  any  LASA  beam  is  the  first 
problem  considered.  This  loss  is  associated  with  a  specified  steering  pattern 
{Uik}  of  preformed  beams.  A  procedure  for  evaluating  such  loss  is  dis¬ 
cussed  below  and  presented  diagrammatic  ally  in  the  flow  chart  of  Figure  84. 


First  a  specified  LASA  beam  j  is  considered.  For  each  subarray  i  the  sub¬ 
array  beam  k  is  determined  so  that  the  vector  magnitude  At  jk  =  -  Ufk  I  iR 

minimal,  and  this  minimal  value  is  designated  by  the  symbol 


The  minimization  is  accomplished  by  calculating  the  set  of  numbers 
{^ijk  (k  =  1.  1  •  •  »  K) }  and  selecting  its  smallest  member  A^*'. 


The  program  now  calculates  the  subarray  loss 
Ak0,  utilizing  Formula  2.31  . 


i  j 


=  A^(f  A^o  )2. 

ij 


corresponding  to 


In  the  above  formula,  f  is  frequency  and  is  a  constant  which  depends  on  the 
subarray  geometry. 
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Input 


{U.},fu(j|,{u(k| 

T 


Current 
Programming 


Sart{  LL.}  In  Decreas- 

ng  order,  where 

LL.  =  Max.{  LL.} 

J  J  J 

- f - 


I 


Future 

Programming  "1 


Calculate  LAS  A 
Lass 

LL. 


Add  1  ta  i 


J 


Output 


A.  List  af  LASA  Beam  Locations  {IL  } 

and  Corresponding  Subarray 
Aiming  Paints{  (T  .  } 

AJ 

B.  List  af  Subarray  Beams  {U^}  Farmed 
far  LASA  Beam  Assignment 

C. l.  List  af  Assigned  Subarray  Losses 

and  Corresponding  Subarrays,  Sub- 
arrays,  Subarray  Beams,  and  LASA 
Array  Lasses 

2.  List  of  Maximum  Subarray  Losses 
and  Corresponding  Subarrays,  Sub- 
array  Beams,  and  LASA  Array  Lasses 

3.  Sorted  List  af  LASA  Losses  {LLj  } 

D.  List  af  Subarray  Beam  Assignments 
and  Assignment  Frequencies 

Far  Each  U  .  Supply  the  Corresponding 
4ka  *a  Delay  Generation 

Program 


.  Yes 

Calculate 

ii 

£ 

5* 

7T 

7T 

Calculate  Sub¬ 
array  Lass 

\j *  kv 


Add  1  ta  i 


ay*> 

/»\ 

✓  LL.  x 

:  jo  > 

Minimal  / 

X  s 


Na 


_ 1 _ 

I  I 

I  f-  ,  I 

j  Reassign  {U^}  | 

I  ' 


Fsgure  84.  Subarray  Assignment  Program 
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There  is  a  loss  at  the  LASA  level,  represented  by  the  symbol  LLj,  which 
occurs  because  of  the  use  of  pre-formed  beams.  This  loss  is  a  function  of  the 

L  subarray  losses:  {X  ...  }  (i  =  1,  .  .  .  ,  L).  Loss  LL.  is  given  by: 

J 


LL.  =  -20  login[T 


/20> 


Subarrays  will  be  identified  by  number  so  that  if  one  or  more  are  deleted,  their 
numbers  will  be  passed  over  in  forming  LLj.  The  purpose  of  such  numbering 
is  to  permit  identification  of  exactly  which  subarrays  are  being  used  in  cases 
where  the  LASA  beam  is  formed  using  fewer  than  the  full  number  of  subarrays. 

The  preceding  calculations  are  iterated  to  obtain  the  sets  of  subarray  losses 
{  ^j£k0  i  }  (j  =  1,  .  .  .  ,  J;  i  =  1,  •  •  .  ,  L)  and  the  set  of  LASA  losses  {LL.} 

(j  =  1.  !  .  .  ,  J)  associated  with  all  LASA  beams  previously  defined.  The  •* 

LASA  losses  {LL-  j  will  be  output  in  decreasing  order.  In  addition,  the  computed 
subarray  losses  ^  will  be  listed.  An  independent  output  of  the  above 

program  will  be  a  list  specifying  the  subarray  beams  to  be  used  in  forming 

each  LASA  beam  U»..  This  list  will  be  supplied  to  the  Time  Delay  Generation 
Program.  ^ 


III. 3. 4  Minimization  of  Loss 

Remaining  problems  that  must  be  formulated  more  precisely  and  for  which 
solutions  are  needed  are  as  follows; 

a.  Formulation  of  criteria  for  acceptable  minimization  of  the  maximum 
loss  LLj0  of  the  set  {LLj  }  (j  =  1,  .  .  .  ,  J). 

b.  Formulation  of  criteria  for  iterated  reassignment  of  the  set  {  Ujfo  }  of 
subarray  assignments.  Such  criteria  should  provide  a  method  of 
feedback  control  maximizing  the  probability  that  iteration  of  the 
outer  loop  of  Figure  84  will  produce  successively  smaller  values 

of  LL;  . 

Jo 
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Appendix  IV 


EVENT  PROCESSING  STUDIES 


Event  processing,  as  outlined  in  Reference  4,  uses  seismometer  data  along 
with  detections  supplied  by  the  Detection  Processor  to  estimate  the  location,  depth, 
magnitude,  and  origin  time  of  an  event.  Inherent  to  the  system  will  be  the  neces¬ 
sary  filtering,  recording,  and  display  processes  with  optional  operator  intervention. 

This  appendix  contains  two  studies  which  support  event  processing.  The  first 
is  concerned  with  an  energy  "density"  method  of  calculating  event  magnitude.  The 
method  is  related  to  the  classical  measure  of  magnitude  for  a  sine  wave  ground 
displacement  model.  The  second  deals  with  the  least  squares  orthogonal  poly¬ 
nomial  fitting  method  used  to  reduce  voluminous  discrete  travel  time  data  to  a 
polynomial  form.  The  polynomials  can  be  used  by  the  Event  Processor  to  deter¬ 
mine  post  P  arrival  windows  and  to  estimate  geographic  locations  of  an  event. 


IV.  1  MAGNITUDE  ESTIMATION 

An  integral  part  of  event  processing  will  be  the  capability7  to  automatically 
estimate  the  parameter  magnitude  which  is  related  to  the  energy  released  at  the 
source  of  the  event.  The  current  classical  technique  uses  manual  measurement 
to  obtain  waveform  amplitude  and  dominant  period.  Computer  extraction  of  these 
parameters  appears  difficult  especially  when  the  signal-to-noise  ratio  is  small. 

In  this  section  an  average  energy  estimation  method  is  presented  and  is 
related  to  the  calculation  of  event  magnitude  for  a  sinusoidal  model  of  ground 
motion.  Furthermore,  the  method  is  expected  to  produce  reasonable  estimates 
of  magnitude  for  non- sinusoidal  waveforms  because,  unlike  the  classical  method, 
all  spectral  components  contribute  to  the  average  energy  estimate. 


IV.  1.2  Classical  Determination  of  Magnitude 

18 

Gutenberg  and  Richter  have  related  the  unified  "magnitude"  of  a  seismic 
event  to  the  compressional  body  wave  amplitude  received  at  a  seismic  station  as: 

10^A 

mb  =  Q  (A>  h)  +  loS10  - Ip-  +  S,  (89) 
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where  Q(A,h)  is  the  depth  distance  function  for  P  wave  using  a  vertical  instru¬ 
ment,  A  is  the  maximum  zero  to  peak  amplitude  of  ground  displacement  in 
the  first  three  cycles  and  measured  in  meters,  T  is  the  dominate  period  in 
seconds,  S  is  the  station  correction,  A  is  the  epicentral  distance  in  degrees 
between  the  event  and  the  seismic  station  and  h  is  the  event  depth  in  kilometers. 

The  above  formula  gives  consistent  results  for  manual  measurements  of  A 
and  T  from  seismograms.  It  is,  however,  not  easily  implemented  on  the  digital 
computer,  nor  does  it  account  for  the  spectral  characteristics  of  the  signal  and 
difficulty  is  expected  when  the  signal  to  noise  ratio  is  low. 


IV.  1.3  Sine  Wave  Ground  Motion  Model  and  Magnitude  Estimate 
19 

Richter  '  notes:  "The  quantity  A/T  has  the  advantage  of  being  simply 
related,  in  theory  at  least,  to  the  kinetic  energy  of  its  wave  train."  Using  this, 
we  assume  that  the  ground  displacement  along  the  axis  of  a  seismometer  result¬ 
ing  from  a  distant  event  is  of  the  form: 

g(t)  =  A  sin  —  t  meters,  (90) 


where  A  is  the  maximum  ground  displacement  in  meters,  T  is  the  period  in 
seconds  and  t  is  the  time  in  seconds.  The  ground  velocity  is: 


v(t)  - 


2  7rA 

T 


cos 


meters /sec. 


(91) 


The  kinetic  energy  of  an  incremental  mass  beneath  the  seismometer  is: 

E’(t)  =  1/2  6  m  v^  (t)  joules,  (92) 

where  6  m  is  the  incremental  mass  in  kilograms.  The  average  kinetic  energy 
normalized  with  respect  to  incremental  mass  is: 

<  E(t)  >  =  TT2  (y)2  joules/kg. 

The  average  energy  for  sinusoidal  ground  displacement  can  now  be  related  to 
event  magnitude. 

To  relate  the  average  energy  measurement  of  the  sinusoidal  model  to 
magnitude,  solve  Equation  (92)  for  |  A  j  and  substitute  in  Equation  (89)  to  give: 

mb  =  0.5  1og1()  |  <  E(t)>  |  +  Q(A,  h)  +  K  +  S,  and  (93) 

K  =  5.502502  =  -log  Q7r  +  6.0  (94) 
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Equation  93  is  valid  for  a  sine  wave  input  and  should  give  a  reasonable  estimate 
of  magnitude  for  other  periodic  waveforms  because  it  includes  the  average  power 
contribution  of  all  spectral  components.  However,  it  does  not  include  the  effects 
of  seismic  noise. 

Because  low  level  seismic  noise  is  always  present,  the  average  energy 
measurement  indicated  in  Equation  (92)  is  a  sum  of  signal  and  noise  energy. 
Assuming  that  the  noise  is  stationary  and  random: 

<Es+n(t)>=  <Es(t)>+  <En(t)>’  (95) 


where  Ec,(t)  is  the  kinetic  energy  of  the  signal  per  unit  mass  and  En(t)  is  the 
kinetic  energy  of  the  noise  per  unit  mass;  so  that 


<E  (t)>  =  <E  ,  (t)> 
s'  ’  s+n'  ' 


<E  (t)> 
n'  ' 


(96) 


Therefore,  average  signal  energy  can  be  found  by  subtracting  the  average  noise 
energy  from  the  measured  average  signal  plus  noise  energy. 

The  magnitude  estimation  process  is  outlined  in  Figure  85.  The  signal  is 
not  truly  periodic,  but  rather  it  is  pulse-like  in  nature  and  has  a  large  amplitude 
for  just  a  few  seconds  of  time.  It  is  during  these  few  seconds  that  the  signal  plus 
noise  average  energy  must  be  estimated.  The  magnitude  estimation  process  is 
accomplished  by  sliding  a  window  along  the  waveform  and  calculating  the  average 
energy  in  the  window.  When  the  window  covers  the  signal,  its  output  is  maximum. 
The  peak  value  is  taken  as  the  estimate  of  the  signal  plus  noise  average  energy. 
Equation  (93)  is  rewritten  as 

mb  =  0.5  1og10  I  <Eg(t)  >  I  p  +  Q(  A,  h)  +  K>  S  (97) 

where  <E^(t)>  is  calculated  from  Equation  (96)  at  the  peak  value  of  <Eg+n(t)>. 

The  signal  plus  noise  energy  estimation  window  should  be  the  same  length  as  the 
significant  signal.  This  is  about  one  or  two  cycles  depending  on  the  signature  of 
the  event.  Similarly,  Equation  (92)  cannot  be  evaluated  exactly  for  seismic  noise 
because  it  is  stationary  for  only  a  few  seconds.  Therefore,  the  average  signal 
energy  and  consequently  the  event  magnitude,  is  based  on  the  two  estimates  of 
the  average  noise  energy  and  the  average  signal  plus  noise  energy. 


IV.  1.4  Sensor  and  Transfer  Function 

The  sensor  system  transfer  function  is  depicted  in  Figure  86.  The  velocity 
transfer  function  is  much  flatter  than  the  displacement  transfer  function  indicating 
that  the  seismometer  is  essentially  a  velocity  measuring  device.  Using  the  instru¬ 
ment's  velocity  measuring  capability,  a  direct  approximation  of  the  quantity  A/T 
can  be  obtained  and  employed  as  indicated  above. 
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Figure  85.  Automatic  Magnitude  Estimation  Process 


(ap)  asuodsay 
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Figure  86.  Instrument  Response 


IV. 2  TRAVEL  TIME  AND  HORIZONTAL  INVERSE  PHASE  VELOCITY 
CHARACTERIZATION 

Seismic  body  wave  travel  time  along  with  the  related  horizontal  inverse 
phase  velocity  are  featured  in  two  areas  of  event  processing. 4  First,  the 
estimation  of  post  P  arrivals  in  both  time  after  the  P  arrival  and  also 
location  in  the  horizontal  inverse  phase  velocity  space  (hereafter  designated 
as  U  space).  The  second  area  of  interest  is  the  conversion  of  U  space  loca¬ 
tions  and  arrival  times  for  P  and  post  P  waves  to  geographic  coordinates  and 
an  estimate  of  event  depth. 

The  purpose  of  this  section  is  to  demonstrate  the  use  of  least  squares 
orthogonal  curve  fitting  in  one  dimension  to  reduce  the  necessary  travel  time 
and  U  space  magnitude  data  to  a  minimum.  Curve  fitting  also  appears  to  auto¬ 
matically  assist  in  interpolation  and  in  obtaining  the  derivative  of  the  discrete 
data. 


IV.2.1  The  Jeffreys-Bullen  Travel  Time  Tables 

1 5 

The  Seismological  Tables  by  Jeffreys  and  Bullen  ,  are  used  to  determine 
event  location  and  time  of  origin  from  event  arrival  times  recorded  at  a  set 
of  seismic  stations.  They  have  been  employed  to  aid  in  identifying  a  method 
for  concise  data  representation  and,  quantitatively,  how  far  such  a  process 
need  be  carried  out. 

The  total  travel  time,  T,  for  a  seismic  body  wave  can  be  found  directly  in 
the  J-B  tables  as  a  two  dimensional  function  of  event  range  and  event  depth. 

The  magnitude  of  the  horizontal  inverse  phase  velocity  vector  in  seconds  per 
kilometer  is  by  definition  the  partial  derivative  of  the  travel  time  with  respect 
to  distance. 

Finding  each  value  of  |u|  from  the  J-B  tables  would  require  calculating 
the  derivative  of  a  low  order  smoothing  function  passed  through  several  imme¬ 
diate  data  points  or  the  physical  storage  of  a  corresponding  table.  Storage  and 
manipulation  of  large  tables  (P  alone  has  over  1400  entries)  can  be  a  burdensome 
process.  Direct  usage  of  the  J-B  tables  to  find  total  travel  time,  T,  and  the 
inverse  horizontal  phase  velocity  from  an  assumed  range  and  depth  is  a  case 
in  point.  While  the  tables  are  consistent  in  having  14  fixed  depths,  the  epicentral 
range  increments  are  not  the  same  for  all  arrival  waves.  In  particular,  the 
change  of  range  from  one  entry  in  the  table  to  the  next  varies  considerably  in  the 
depth  allowance  tables. 

As  an  aid  in  identifying  later  arrivals  of  an  event,  there  is  a  need  to  compute 
theoretical  travel  time  differences  between  an  assumed  P  wave  and  the  later 
arrivals,  and  the  utilization  of  approximating  analytical  functions  must  be 
considered. 
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IV. 2. 2  Table  Approximation  by  Classical  Least  Squares 

Since  direct  usage  of  the  J-B  tables  would  be  computationally  unwieldy  in 
event  processing,  attention  was  given  to  the  classical  curve  fitting  of  data  points 
in  the  least  squares  sense.  The  fit  to  P  wave  data  of  surface  depth  deteriorated 
beyond  fourth  order  despite  the  use  of  double  precision  (16  digit)  arithmetic 
computer  operations.  This  was  not  surprising,  since  the  classical  method 
requires  the  inversion  of  a  matrix  and  involves  the  notorious  principal  minor  of 
the  well  known  Hilbert  matrix^. 

Curve  fitting,  by  the  use  of  polynomials  which  are  orthogonal  under  summa¬ 
tion,  was  deemed  the  next  logical  step  since  their  use  does  not  require  the  inver¬ 
sion  of  a  matrix.  While  tabulations  of  precalculated  orthogonal  polynomials  are 
available*^,  their  application  requires  equal  spacing  of  the  independent  variable 
increments. 

Since  the  J-B  tables  have  an  unequally  spaced  variable  of  epicentral  range, 
Forsythe's  regression  form  of  orthogonal  polynomial  curve  fitting  was  adopted“^*>22. 
A  brief  description  of  this  method  follows. 

Given  data  for  (m)  points  (xk,  yk),  k  =  1,  2,  .  .  .  ,  m;  express  (y)  as  a  func¬ 
tion  of  (x)  by  the  use  of  orthogonal  polynomials  of  degree  0  through  degree  (n), 
where  (n  <  m  -  1).  Thus,  find  bg,  bj,  b2»  .....  bn  so  that: 

y(x)  =  t>0P0(x)  +  bjP^x)  +  .  .  .  +  bnPn(x),  (98) 

where  P^(x)  is  a  polynomial  of  degree  (i)  and  y(x)  best  fits  the  data  in  the  sense 
of  least  squares: 


-  y  (xk)l 


2  . 

is  a  minimum, 


where  yk  is  the  given  tabular  value  and  y(xk)  is  the  corresponding  calculated  value. 


After  y(x)  has  been  expressed  in  the  form  98,  it  can  be  converted  to  a  poly¬ 
nomial  in  powers  of  x  since  each  Pi(x)  is  a  polynomial  in  x: 


2  n 

y(x)  =  an  +  a  x  +  a  x  +  .  .  .  +  a  x 


(99) 
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The  coefficients  b((,  b^,  b.,,  .  .  .  ,  bn  may  be  found  from  the  formula: 


b.  = 

1 


m 


k=i 


(100) 


The  orthogonal  polynomials  are  generated  by  a  recursive  formula  as 
follows: 


P^x)  =  0 


P0(x) 


P.+1(x)  =  X  P.(x)  -  fci+1P.(x)  -  ^  Pj.^x). 


(101) 


The  coefficients  and  (3 .  are  chosen  so  that  the  polynomials  P^,  P^,  P 
,  Pn  are  orthogonal  under  summation: 

m 

\  pi(xk)pj(xk)  =  0  for  a11  (i  ^  J)-  <102) 

k-1 


The  formulas  for  a.  1  and  B.  are  as  follows: 

l+l  '  l 


m 


I  xk(Pi 


(Xu))‘ 


a.  ,  = 


k=l 


1+1  m 


X*. 

k-1 


(xi ,)) 


(103) 
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m 


(pi(xk>r 


p  .  = 


k=l 


i  m 


(pi_l(xk))' 


<=1 


(104) 


Forsythe  suggests  transforming  the  independent  variable  x  to  the  range  -2. 
to  +2.  to  avoid  possible  overflow.  This  transformation  also  permits  fewer 
digits  to  be  used  for  each  coefficient  a^  for  a  specific  accuracy  of  fit. 

IV. 2. 3  Results  of  Least  Squares  Orthogonal  Polynomial  Fit  to  the 
Jeffreys-Bullen  Tables 

The  Forsythe  regression  form  of  orthogonal  polynominal  curve  fitting  was 
used  to  reduce  the  tables  to  manageable  polynomials.  The  errors  of  fitting  the 
travel  time,  T,  for  the  P  wave,  PcP  wave,  and  pP-P  wave  difference  arc  evident 
on  figures  87  to  90. 

Figure  87  is  primarily  a  plot  of  the  error  of  fit  in  seconds  versus  the  order 
of  the  polynomial  fit.  Two  indications  of  fitting  error  are  given:  the  standard 
deviation  (RMS  error)  and  the  maximum  deviation.  Both  are  plotted  as  a  range 
covering  the  14  depths.  Using  |  U  |  -  0.04  sec /Km  at  30°  distance  as  a  worst  case, 
a  timing  error  of  one  second  is  equivalent  to  a  25 Km  error  in  location  and  is 
called  the  ’location  error  bound”  on  the  graph.  Q  is  the  quantization  level  of  the 
original  J-B  data  and  a  is  the  expected  RMS  error  in  the  original  data  caused 
by  quantization. 

The  ninth  order  fit  for  T  seems  best  because  the  maximum  location  error  is 
only  3  Km  and  the  RMS  error  value  of  the  fit  was  bounded  by  a. 

Figure  88  is  very  similar  to  Figure  87  except  that  the  wave  in  question  is 
PcP  and  for  an  event  at  the  surface  only.  Using  |u|  =  0.023  sec/Km  at  30°  dis¬ 
tance  as  a  worst  case,  a  timing  error  of  one  second  is  equivalent  to  a  43.5  Km 
error  in  location.  An  eighth  order  fit  looks  best  for  T. 

Figure  89  shows  the  error  of  fit  for  the  13  depth  correction  polynomials  for 
PcP.  The  d/RQ  is  a  standard  form  for  showing  depth  in  a  compact  form  and  can 
be  viewed  as  ratio  of  the  event  depth  below  the  crust  (the  crust  is  about  33  Km 
thick)  to  the  radius  of  the  earth  minus  the  crustal  layer.  Thus,  a  depth  of 
33  Km,  as  measured  from  the  surface  of  the  earth,  becomes  equal  to  zero  in 
d/Rg  format.  The  worst  case  for  an  error  in  the  estimation  of  depth  due  to  a 
timing  error  occurs  at  80°  range  and  d/Rp  =  0.11.  Here,  only  5.1  seconds 
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Order  of  Polynomial 

Moximum  Deviotion  Over  All  Depths  (0.  to 794  Km) 

Standord  Deviotion  Over  All  depths  (0.  to 794  Km) 

Is  Data  Quantization  Level  of  Jeffreys-Bul  len  P  Table 

Is  Expected  RMS  Error  in  Dato  Caused  By  Quantization  1 


Figure  8 70  Least  Squares  Fit  for  P  Wave  at  14  Depths 
for  25°  -  105°  Range 
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Locotion  Error  Bound  in 


Error  of  Fit  in  sec 


Order  of  Polynomial 


^  Maximum  Deviation 
O  Standard  Deviation 

Q  Is  Data  Quantization  Level  of  Jeffreys-Bul len  P^P  Table 
O  Is  Expected  RMS  Error  in  Data  Caused  by  Quantization 


Figure  880  Least  Squares  for  PcP  Wave  at  Surface  Depths 
for  25°  *  100°  Range 
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Location  Error  Bound  in  Km 
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Figure  89„  Leas*  Squares  Fbt  for  PcP  Wave  at  13  Depths 
for  0°  -  80°  Range 
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Fugure  90o  Least  Squares  F!t  for  pP-P  at  13  Depths 
for  30°  -  100°  Range 
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Depth  Error  Bound  in 


separates  the  two  depth  columns  (0.11  and  0.12)  in  the  J-B  tables.  This  results 
in  an  upper  bound  of  12.4  km  error  in  depth  for  a  one  second  timing  error. 

Figure  90  is  similar  in  design^to  Figure  89  and  shows  the  error  of  fit  for 
the  pP-P  time  difference.  Since  |  U  |  is  the  derivative  of  travel  time  with  respect 
to  distance,  the  order  of  the  |u  |  polynomial  will  be  one  order  less  than  the  proper 
fit  for  T, 


IV. 2.4  Conclusions  and  Recommendations 

Forsythe's  regression  form  of  orthogonal  polynomial  curve  fitting  outper¬ 
forms  the  classical  least  squares  method  and  appears  to  be  quite  useful  in  fitting 
one  dimensional  travel  time  data.  While  taking  the  derivative  of  the  polynomial  fit 
to  obtain  |  0]  as  a  function  of  range  seems  better  than  the  derivative  of  discrete 
data,  it  may  offer  only  marginal  accuracy  in  geographic  conversion  yet  be  per¬ 
fectly  acceptable  for  post  P  detection  windows  in  U  space.  Fitting  T  and  |  u|  data 
obtained  at  LASA  should  be  no  problem  since  they  are  direct  measurements.  It  is 
recommended  that  a  two  dimensional  surface  fitting  method  be  developed  so  that 
both  range  and  depth  are  independent  parameters. 
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